libsim Versione 7.2.1
|
◆ georef_coord_array_compute_bbox()
Compute the bounding box of each shape in georef_coord_array object. The bounding box is computed and stored in the object, it is used by the inside() function for speedup; after it is computed the object cannot be changed, otherwise the bounding box will not be valid.
Definizione alla linea 1035 del file georef_coord_class.F90. 1036! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1037! authors:
1038! Davide Cesari <dcesari@arpa.emr.it>
1039! Paolo Patruno <ppatruno@arpa.emr.it>
1040
1041! This program is free software; you can redistribute it and/or
1042! modify it under the terms of the GNU General Public License as
1043! published by the Free Software Foundation; either version 2 of
1044! the License, or (at your option) any later version.
1045
1046! This program is distributed in the hope that it will be useful,
1047! but WITHOUT ANY WARRANTY; without even the implied warranty of
1048! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1049! GNU General Public License for more details.
1050
1051! You should have received a copy of the GNU General Public License
1052! along with this program. If not, see <http://www.gnu.org/licenses/>.
1053#include "config.h"
1054
1072USE geo_proj_class
1073#ifdef HAVE_SHAPELIB
1074USE shapelib
1075#endif
1076IMPLICIT NONE
1077
1083 PRIVATE
1084 DOUBLE PRECISION :: x=dmiss, y=dmiss
1086
1089
1095 PRIVATE
1096 INTEGER,ALLOCATABLE :: parts(:)
1097 TYPE(georef_coord),ALLOCATABLE :: coord(:)
1098 INTEGER :: topo=imiss
1099 TYPE(geo_proj) :: proj
1100 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
1101 LOGICAL :: bbox_updated=.false.
1103
1104INTEGER,PARAMETER :: georef_coord_array_point = 1
1105INTEGER,PARAMETER :: georef_coord_array_arc = 3
1106INTEGER,PARAMETER :: georef_coord_array_polygon = 5
1107INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
1108
1109
1114 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1115END INTERFACE
1116
1119 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1120END INTERFACE
1121
1124 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
1125END INTERFACE
1126
1127INTERFACE compute_bbox
1128 MODULE PROCEDURE georef_coord_array_compute_bbox
1129END INTERFACE
1130
1132INTERFACE OPERATOR (==)
1133 MODULE PROCEDURE georef_coord_eq
1134END INTERFACE
1135
1137INTERFACE OPERATOR (/=)
1138 MODULE PROCEDURE georef_coord_ne
1139END INTERFACE
1140
1143INTERFACE OPERATOR (>=)
1144 MODULE PROCEDURE georef_coord_ge
1145END INTERFACE
1146
1149INTERFACE OPERATOR (<=)
1150 MODULE PROCEDURE georef_coord_le
1151END INTERFACE
1152
1153#ifdef HAVE_SHAPELIB
1154
1156INTERFACE import
1157 MODULE PROCEDURE arrayof_georef_coord_array_import
1158END INTERFACE
1159
1163 MODULE PROCEDURE arrayof_georef_coord_array_export
1164END INTERFACE
1165#endif
1166
1170 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1171END INTERFACE
1172
1176 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1177END INTERFACE
1178
1181 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1182END INTERFACE
1183
1186 MODULE PROCEDURE georef_coord_dist
1187END INTERFACE
1188
1189#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1190#define ARRAYOF_TYPE arrayof_georef_coord_array
1191!define ARRAYOF_ORIGEQ 0
1192#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1193#include "arrayof_pre.F90"
1194! from arrayof
1196
1197PRIVATE
1199 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1200 georef_coord_array_polygon, georef_coord_array_multipoint, &
1202#ifdef HAVE_SHAPELIB
1204#endif
1206 georef_coord_new, georef_coord_array_new
1207
1208CONTAINS
1209
1210#include "arrayof_post.F90"
1211
1212! ===================
1213! == georef_coord ==
1214! ===================
1218FUNCTION georef_coord_new(x, y) RESULT(this)
1219DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1220DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1221TYPE(georef_coord) :: this
1222
1225
1226END FUNCTION georef_coord_new
1227
1228
1229SUBROUTINE georef_coord_delete(this)
1230TYPE(georef_coord),INTENT(inout) :: this
1231
1232this%x = dmiss
1233this%y = dmiss
1234
1235END SUBROUTINE georef_coord_delete
1236
1237
1238ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1239TYPE(georef_coord),INTENT(in) :: this
1240LOGICAL :: res
1241
1242res = .NOT. this == georef_coord_miss
1243
1244END FUNCTION georef_coord_c_e
1245
1246
1253ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1254TYPE(georef_coord),INTENT(in) :: this
1255DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1256DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1257
1258IF (PRESENT(x)) x = this%x
1259IF (PRESENT(y)) y = this%y
1260
1261END SUBROUTINE georef_coord_getval
1262
1263
1272ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1273TYPE(georef_coord),INTENT(in) :: this
1274TYPE(geo_proj),INTENT(in) :: proj
1275DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1276DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1277DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1278DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1279
1280DOUBLE PRECISION :: llon, llat
1281
1282IF (PRESENT(x)) x = this%x
1283IF (PRESENT(y)) y = this%y
1284IF (PRESENT(lon) .OR. present(lat)) THEN
1286 IF (PRESENT(lon)) lon = llon
1287 IF (PRESENT(lat)) lat = llat
1288ENDIF
1289
1290END SUBROUTINE georef_coord_proj_getval
1291
1292
1293! document and improve
1294ELEMENTAL FUNCTION getlat(this)
1295TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1296DOUBLE PRECISION :: getlat ! latitudine geografica
1297
1298getlat = this%y ! change!!!
1299
1300END FUNCTION getlat
1301
1302! document and improve
1303ELEMENTAL FUNCTION getlon(this)
1304TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1305DOUBLE PRECISION :: getlon ! longitudine geografica
1306
1307getlon = this%x ! change!!!
1308
1309END FUNCTION getlon
1310
1311
1312ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1313TYPE(georef_coord),INTENT(IN) :: this, that
1314LOGICAL :: res
1315
1316res = (this%x == that%x .AND. this%y == that%y)
1317
1318END FUNCTION georef_coord_eq
1319
1320
1321ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1322TYPE(georef_coord),INTENT(IN) :: this, that
1323LOGICAL :: res
1324
1325res = (this%x >= that%x .AND. this%y >= that%y)
1326
1327END FUNCTION georef_coord_ge
1328
1329
1330ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1331TYPE(georef_coord),INTENT(IN) :: this, that
1332LOGICAL :: res
1333
1334res = (this%x <= that%x .AND. this%y <= that%y)
1335
1336END FUNCTION georef_coord_le
1337
1338
1339ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1340TYPE(georef_coord),INTENT(IN) :: this, that
1341LOGICAL :: res
1342
1343res = .NOT.(this == that)
1344
1345END FUNCTION georef_coord_ne
1346
1347
1353SUBROUTINE georef_coord_read_unit(this, unit)
1354TYPE(georef_coord),INTENT(out) :: this
1355INTEGER, INTENT(in) :: unit
1356
1357CALL georef_coord_vect_read_unit((/this/), unit)
1358
1359END SUBROUTINE georef_coord_read_unit
1360
1361
1367SUBROUTINE georef_coord_vect_read_unit(this, unit)
1368TYPE(georef_coord) :: this(:)
1369INTEGER, INTENT(in) :: unit
1370
1371CHARACTER(len=40) :: form
1372INTEGER :: i
1373
1374INQUIRE(unit, form=form)
1375IF (form == 'FORMATTED') THEN
1376 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1377!TODO bug gfortran compiler !
1378!missing values are unredeable when formatted
1379ELSE
1380 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1381ENDIF
1382
1383END SUBROUTINE georef_coord_vect_read_unit
1384
1385
1390SUBROUTINE georef_coord_write_unit(this, unit)
1391TYPE(georef_coord),INTENT(in) :: this
1392INTEGER, INTENT(in) :: unit
1393
1394CALL georef_coord_vect_write_unit((/this/), unit)
1395
1396END SUBROUTINE georef_coord_write_unit
1397
1398
1403SUBROUTINE georef_coord_vect_write_unit(this, unit)
1404TYPE(georef_coord),INTENT(in) :: this(:)
1405INTEGER, INTENT(in) :: unit
1406
1407CHARACTER(len=40) :: form
1408INTEGER :: i
1409
1410INQUIRE(unit, form=form)
1411IF (form == 'FORMATTED') THEN
1412 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1413!TODO bug gfortran compiler !
1414!missing values are unredeable when formatted
1415ELSE
1416 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1417ENDIF
1418
1419END SUBROUTINE georef_coord_vect_write_unit
1420
1421
1424FUNCTION georef_coord_dist(this, that) RESULT(dist)
1426TYPE(georef_coord), INTENT (IN) :: this
1427TYPE(georef_coord), INTENT (IN) :: that
1428DOUBLE PRECISION :: dist
1429
1430DOUBLE PRECISION :: x,y
1431! Distanza approssimata, valida per piccoli angoli
1432
1433x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1434y = (this%y-that%y)
1435dist = sqrt(x**2 + y**2)*degrad*rearth
1436
1437END FUNCTION georef_coord_dist
1438
1439
1445FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1446TYPE(georef_coord),INTENT(IN) :: this
1447TYPE(georef_coord),INTENT(IN) :: coordmin
1448TYPE(georef_coord),INTENT(IN) :: coordmax
1449LOGICAL :: res
1450
1451res = (this >= coordmin .AND. this <= coordmax)
1452
1453END FUNCTION georef_coord_inside_rectang
1454
1455
1456! ========================
1457! == georef_coord_array ==
1458! ========================
1464FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1465DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1466DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1467INTEGER,INTENT(in),OPTIONAL :: topo
1468TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1469TYPE(georef_coord_array) :: this
1470
1471INTEGER :: lsize
1472
1473IF (PRESENT(x) .AND. PRESENT(y)) THEN
1474 lsize = min(SIZE(x), SIZE(y))
1475 ALLOCATE(this%coord(lsize))
1476 this%coord(1:lsize)%x = x(1:lsize)
1477 this%coord(1:lsize)%y = y(1:lsize)
1478ENDIF
1479this%topo = optio_l(topo)
1481
1482END FUNCTION georef_coord_array_new
1483
1484
1485SUBROUTINE georef_coord_array_delete(this)
1486TYPE(georef_coord_array),INTENT(inout) :: this
1487
1488TYPE(georef_coord_array) :: lobj
1489
1490this = lobj
1491
1492END SUBROUTINE georef_coord_array_delete
1493
1494
1495ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1496TYPE(georef_coord_array),INTENT(in) :: this
1497LOGICAL :: res
1498
1499res = ALLOCATED(this%coord)
1500
1501END FUNCTION georef_coord_array_c_e
1502
1503
1508SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1509TYPE(georef_coord_array),INTENT(in) :: this
1510DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1511DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1512! allocatable per vedere di nascosto l'effetto che fa
1513INTEGER,OPTIONAL,INTENT(out) :: topo
1514TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1515
1516
1517IF (PRESENT(x)) THEN
1518 IF (ALLOCATED(this%coord)) THEN
1519 x = this%coord%x
1520 ENDIF
1521ENDIF
1522IF (PRESENT(y)) THEN
1523 IF (ALLOCATED(this%coord)) THEN
1524 y = this%coord%y
1525 ENDIF
1526ENDIF
1527IF (PRESENT(topo)) topo = this%topo
1529
1530END SUBROUTINE georef_coord_array_getval
1531
1532
1538SUBROUTINE georef_coord_array_compute_bbox(this)
1539TYPE(georef_coord_array),INTENT(inout) :: this
1540
1541IF (ALLOCATED(this%coord)) THEN
1542 this%bbox(1)%x = minval(this%coord(:)%x)
1543 this%bbox(1)%y = minval(this%coord(:)%y)
1544 this%bbox(2)%x = maxval(this%coord(:)%x)
1545 this%bbox(2)%y = maxval(this%coord(:)%y)
1546 this%bbox_updated = .true.
1547ENDIF
1548
1549END SUBROUTINE georef_coord_array_compute_bbox
1550
1551#ifdef HAVE_SHAPELIB
1552! internal method for importing a single shape
1553SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1554TYPE(georef_coord_array),INTENT(OUT) :: this
1555TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1556INTEGER,INTENT(IN) :: nshp
1557
1558TYPE(shpobject) :: shpobj
1559
1560! read shape object
1561shpobj = shpreadobject(shphandle, nshp)
1562IF (.NOT.shpisnull(shpobj)) THEN
1563! import it in georef_coord object
1564 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1565 topo=shpobj%nshptype)
1566 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1567 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1568 ELSE IF (ALLOCATED(this%parts)) THEN
1569 DEALLOCATE(this%parts)
1570 ENDIF
1571 CALL shpdestroyobject(shpobj)
1572 CALL compute_bbox(this)
1573ENDIF
1574
1575
1576END SUBROUTINE georef_coord_array_import
1577
1578
1579! internal method for exporting a single shape
1580SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1581TYPE(georef_coord_array),INTENT(in) :: this
1582TYPE(shpfileobject),INTENT(inout) :: shphandle
1583INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1584
1585INTEGER :: i
1586TYPE(shpobject) :: shpobj
1587
1588IF (ALLOCATED(this%coord)) THEN
1589 IF (ALLOCATED(this%parts)) THEN
1590 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1591 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1592 ELSE
1593 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1594 this%coord(:)%x, this%coord(:)%y)
1595 ENDIF
1596ELSE
1597 RETURN
1598ENDIF
1599
1600IF (.NOT.shpisnull(shpobj)) THEN
1601 i = shpwriteobject(shphandle, nshp, shpobj)
1602 CALL shpdestroyobject(shpobj)
1603ENDIF
1604
1605END SUBROUTINE georef_coord_array_export
1606
1607
1618SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1619TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1620CHARACTER(len=*),INTENT(in) :: shpfile
1621
1622REAL(kind=fp_d) :: minb(4), maxb(4)
1623INTEGER :: i, ns, shptype, dbfnf, dbfnr
1624TYPE(shpfileobject) :: shphandle
1625
1626shphandle = shpopen(trim(shpfile), 'rb')
1627IF (shpfileisnull(shphandle)) THEN
1628 ! log here
1629 CALL raise_error()
1630 RETURN
1631ENDIF
1632
1633! get info about file
1634CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1635IF (ns > 0) THEN ! allocate and read the object
1637 DO i = 1, ns
1638 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1639 ENDDO
1640ENDIF
1641
1642CALL shpclose(shphandle)
1643! pack object to save memory
1645
1646END SUBROUTINE arrayof_georef_coord_array_import
1647
1648
1654SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1655TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1656CHARACTER(len=*),INTENT(in) :: shpfile
1657
1658INTEGER :: i
1659TYPE(shpfileobject) :: shphandle
1660
1661IF (this%arraysize > 0) THEN
1662 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1663ELSE
1664 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1665ENDIF
1666IF (shpfileisnull(shphandle)) THEN
1667 ! log here
1668 CALL raise_error()
1669 RETURN
1670ENDIF
1671
1672DO i = 1, this%arraysize
1673 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1674ENDDO
1675
1676CALL shpclose(shphandle)
1677
1678END SUBROUTINE arrayof_georef_coord_array_export
1679#endif
1680
1692FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1693TYPE(georef_coord), INTENT(IN) :: this
1694TYPE(georef_coord_array), INTENT(IN) :: poly
1695LOGICAL :: inside
1696
1697INTEGER :: i
1698
1699inside = .false.
1701IF (.NOT.ALLOCATED(poly%coord)) RETURN
1702! if outside bounding box stop here
1703IF (poly%bbox_updated) THEN
1704 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1705ENDIF
1706
1707IF (ALLOCATED(poly%parts)) THEN
1708 DO i = 1, SIZE(poly%parts)-1
1710 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1711 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1712 ENDDO
1713 IF (SIZE(poly%parts) > 0) THEN ! safety check
1715 poly%coord(poly%parts(i)+1:)%x, &
1716 poly%coord(poly%parts(i)+1:)%y)
1717 ENDIF
1718
1719ELSE
1720 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1721 inside = pointinpoly(this%x, this%y, &
1722 poly%coord(:)%x, poly%coord(:)%y)
1723ENDIF
1724
1725CONTAINS
1726
1727FUNCTION pointinpoly(x, y, px, py)
1728DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1729LOGICAL :: pointinpoly
1730
1731INTEGER :: i, j, starti
1732
1733pointinpoly = .false.
1734
1735IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1736 starti = 2
1737 j = 1
1738ELSE ! unclosed polygon
1739 starti = 1
1740 j = SIZE(px)
1741ENDIF
1742DO i = starti, SIZE(px)
1743 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1744 (py(j) <= y .AND. y < py(i))) THEN
1745 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1746 pointinpoly = .NOT. pointinpoly
1747 ENDIF
1748 ENDIF
1749 j = i
1750ENDDO
1751
1752END FUNCTION pointinpoly
1753
1754END FUNCTION georef_coord_inside
1755
1756
1757
Compute forward coordinate transformation from geographical system to projected system. Definition geo_proj_class.F90:289 Compute backward coordinate transformation from projected system to geographical system. Definition geo_proj_class.F90:295 Quick method to append an element to the array. Definition georef_coord_class.F90:395 Compute the distance in m between two points. Definition georef_coord_class.F90:338 Export an array of georef_coord_array objects to a file in ESRI/Shapefile format. Definition georef_coord_class.F90:315 Methods for returning the value of object members. Definition georef_coord_class.F90:276 Import an array of georef_coord_array objects from a file in ESRI/Shapefile format. Definition georef_coord_class.F90:309 Method for inserting elements of the array at a desired position. Definition georef_coord_class.F90:386 Determine whether a point lies inside a polygon or a rectangle. Definition georef_coord_class.F90:333 Method for packing the array object reducing at a minimum the memory occupation, without destroying i... Definition georef_coord_class.F90:418 Read a single georef_coord object or an array of georef_coord objects from a Fortran FORMATTED or UNF... Definition georef_coord_class.F90:322 Method for removing elements of the array at a desired position. Definition georef_coord_class.F90:401 Write a single georef_coord object or an array of georef_coord objects to a Fortran FORMATTED or UNFO... Definition georef_coord_class.F90:328 Generic subroutine for checking OPTIONAL parameters. Definition optional_values.f90:36 This module defines objects describing georeferenced sparse points possibly with topology and project... Definition georef_coord_class.F90:221 Definitions of constants and functions for working with missing values. Definition missing_values.f90:50 Module for quickly interpreting the OPTIONAL parameters passed to a subprogram. Definition optional_values.f90:28 Derived type defining a one-dimensional array of georeferenced points with an associated topology (is... Definition georef_coord_class.F90:247 Derive type defining a single georeferenced point, either in geodetic or in projected coordinates. Definition georef_coord_class.F90:235 |