libsim Versione 7.1.11
|
◆ 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 1041 del file georef_coord_class.F90. 1042! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1043! authors:
1044! Davide Cesari <dcesari@arpa.emr.it>
1045! Paolo Patruno <ppatruno@arpa.emr.it>
1046
1047! This program is free software; you can redistribute it and/or
1048! modify it under the terms of the GNU General Public License as
1049! published by the Free Software Foundation; either version 2 of
1050! the License, or (at your option) any later version.
1051
1052! This program is distributed in the hope that it will be useful,
1053! but WITHOUT ANY WARRANTY; without even the implied warranty of
1054! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1055! GNU General Public License for more details.
1056
1057! You should have received a copy of the GNU General Public License
1058! along with this program. If not, see <http://www.gnu.org/licenses/>.
1059#include "config.h"
1060
1078USE geo_proj_class
1079#ifdef HAVE_SHAPELIB
1080USE shapelib
1081#endif
1082IMPLICIT NONE
1083
1089 PRIVATE
1090 DOUBLE PRECISION :: x=dmiss, y=dmiss
1092
1095
1101 PRIVATE
1102 INTEGER,ALLOCATABLE :: parts(:)
1103 TYPE(georef_coord),ALLOCATABLE :: coord(:)
1104 INTEGER :: topo=imiss
1105 TYPE(geo_proj) :: proj
1106 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
1107 LOGICAL :: bbox_updated=.false.
1109
1110INTEGER,PARAMETER :: georef_coord_array_point = 1
1111INTEGER,PARAMETER :: georef_coord_array_arc = 3
1112INTEGER,PARAMETER :: georef_coord_array_polygon = 5
1113INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
1114
1115
1120 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1121END INTERFACE
1122
1125 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1126END INTERFACE
1127
1130 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
1131END INTERFACE
1132
1133INTERFACE compute_bbox
1134 MODULE PROCEDURE georef_coord_array_compute_bbox
1135END INTERFACE
1136
1138INTERFACE OPERATOR (==)
1139 MODULE PROCEDURE georef_coord_eq
1140END INTERFACE
1141
1143INTERFACE OPERATOR (/=)
1144 MODULE PROCEDURE georef_coord_ne
1145END INTERFACE
1146
1149INTERFACE OPERATOR (>=)
1150 MODULE PROCEDURE georef_coord_ge
1151END INTERFACE
1152
1155INTERFACE OPERATOR (<=)
1156 MODULE PROCEDURE georef_coord_le
1157END INTERFACE
1158
1159#ifdef HAVE_SHAPELIB
1160
1162INTERFACE import
1163 MODULE PROCEDURE arrayof_georef_coord_array_import
1164END INTERFACE
1165
1169 MODULE PROCEDURE arrayof_georef_coord_array_export
1170END INTERFACE
1171#endif
1172
1176 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1177END INTERFACE
1178
1182 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1183END INTERFACE
1184
1187 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1188END INTERFACE
1189
1192 MODULE PROCEDURE georef_coord_dist
1193END INTERFACE
1194
1195#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1196#define ARRAYOF_TYPE arrayof_georef_coord_array
1197!define ARRAYOF_ORIGEQ 0
1198#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1199#include "arrayof_pre.F90"
1200! from arrayof
1202
1203PRIVATE
1205 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1206 georef_coord_array_polygon, georef_coord_array_multipoint, &
1208#ifdef HAVE_SHAPELIB
1210#endif
1212 georef_coord_new, georef_coord_array_new
1213
1214CONTAINS
1215
1216#include "arrayof_post.F90"
1217
1218! ===================
1219! == georef_coord ==
1220! ===================
1224FUNCTION georef_coord_new(x, y) RESULT(this)
1225DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1226DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1227TYPE(georef_coord) :: this
1228
1231
1232END FUNCTION georef_coord_new
1233
1234
1235SUBROUTINE georef_coord_delete(this)
1236TYPE(georef_coord),INTENT(inout) :: this
1237
1238this%x = dmiss
1239this%y = dmiss
1240
1241END SUBROUTINE georef_coord_delete
1242
1243
1244ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1245TYPE(georef_coord),INTENT(in) :: this
1246LOGICAL :: res
1247
1248res = .NOT. this == georef_coord_miss
1249
1250END FUNCTION georef_coord_c_e
1251
1252
1259ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1260TYPE(georef_coord),INTENT(in) :: this
1261DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1262DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1263
1264IF (PRESENT(x)) x = this%x
1265IF (PRESENT(y)) y = this%y
1266
1267END SUBROUTINE georef_coord_getval
1268
1269
1278ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1279TYPE(georef_coord),INTENT(in) :: this
1280TYPE(geo_proj),INTENT(in) :: proj
1281DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1282DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1283DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1284DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1285
1286DOUBLE PRECISION :: llon, llat
1287
1288IF (PRESENT(x)) x = this%x
1289IF (PRESENT(y)) y = this%y
1290IF (PRESENT(lon) .OR. present(lat)) THEN
1292 IF (PRESENT(lon)) lon = llon
1293 IF (PRESENT(lat)) lat = llat
1294ENDIF
1295
1296END SUBROUTINE georef_coord_proj_getval
1297
1298
1299! document and improve
1300ELEMENTAL FUNCTION getlat(this)
1301TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1302DOUBLE PRECISION :: getlat ! latitudine geografica
1303
1304getlat = this%y ! change!!!
1305
1306END FUNCTION getlat
1307
1308! document and improve
1309ELEMENTAL FUNCTION getlon(this)
1310TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1311DOUBLE PRECISION :: getlon ! longitudine geografica
1312
1313getlon = this%x ! change!!!
1314
1315END FUNCTION getlon
1316
1317
1318ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1319TYPE(georef_coord),INTENT(IN) :: this, that
1320LOGICAL :: res
1321
1322res = (this%x == that%x .AND. this%y == that%y)
1323
1324END FUNCTION georef_coord_eq
1325
1326
1327ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1328TYPE(georef_coord),INTENT(IN) :: this, that
1329LOGICAL :: res
1330
1331res = (this%x >= that%x .AND. this%y >= that%y)
1332
1333END FUNCTION georef_coord_ge
1334
1335
1336ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1337TYPE(georef_coord),INTENT(IN) :: this, that
1338LOGICAL :: res
1339
1340res = (this%x <= that%x .AND. this%y <= that%y)
1341
1342END FUNCTION georef_coord_le
1343
1344
1345ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1346TYPE(georef_coord),INTENT(IN) :: this, that
1347LOGICAL :: res
1348
1349res = .NOT.(this == that)
1350
1351END FUNCTION georef_coord_ne
1352
1353
1359SUBROUTINE georef_coord_read_unit(this, unit)
1360TYPE(georef_coord),INTENT(out) :: this
1361INTEGER, INTENT(in) :: unit
1362
1363CALL georef_coord_vect_read_unit((/this/), unit)
1364
1365END SUBROUTINE georef_coord_read_unit
1366
1367
1373SUBROUTINE georef_coord_vect_read_unit(this, unit)
1374TYPE(georef_coord) :: this(:)
1375INTEGER, INTENT(in) :: unit
1376
1377CHARACTER(len=40) :: form
1378INTEGER :: i
1379
1380INQUIRE(unit, form=form)
1381IF (form == 'FORMATTED') THEN
1382 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1383!TODO bug gfortran compiler !
1384!missing values are unredeable when formatted
1385ELSE
1386 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1387ENDIF
1388
1389END SUBROUTINE georef_coord_vect_read_unit
1390
1391
1396SUBROUTINE georef_coord_write_unit(this, unit)
1397TYPE(georef_coord),INTENT(in) :: this
1398INTEGER, INTENT(in) :: unit
1399
1400CALL georef_coord_vect_write_unit((/this/), unit)
1401
1402END SUBROUTINE georef_coord_write_unit
1403
1404
1409SUBROUTINE georef_coord_vect_write_unit(this, unit)
1410TYPE(georef_coord),INTENT(in) :: this(:)
1411INTEGER, INTENT(in) :: unit
1412
1413CHARACTER(len=40) :: form
1414INTEGER :: i
1415
1416INQUIRE(unit, form=form)
1417IF (form == 'FORMATTED') THEN
1418 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1419!TODO bug gfortran compiler !
1420!missing values are unredeable when formatted
1421ELSE
1422 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1423ENDIF
1424
1425END SUBROUTINE georef_coord_vect_write_unit
1426
1427
1430FUNCTION georef_coord_dist(this, that) RESULT(dist)
1432TYPE(georef_coord), INTENT (IN) :: this
1433TYPE(georef_coord), INTENT (IN) :: that
1434DOUBLE PRECISION :: dist
1435
1436DOUBLE PRECISION :: x,y
1437! Distanza approssimata, valida per piccoli angoli
1438
1439x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1440y = (this%y-that%y)
1441dist = sqrt(x**2 + y**2)*degrad*rearth
1442
1443END FUNCTION georef_coord_dist
1444
1445
1451FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1452TYPE(georef_coord),INTENT(IN) :: this
1453TYPE(georef_coord),INTENT(IN) :: coordmin
1454TYPE(georef_coord),INTENT(IN) :: coordmax
1455LOGICAL :: res
1456
1457res = (this >= coordmin .AND. this <= coordmax)
1458
1459END FUNCTION georef_coord_inside_rectang
1460
1461
1462! ========================
1463! == georef_coord_array ==
1464! ========================
1470FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1471DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1472DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1473INTEGER,INTENT(in),OPTIONAL :: topo
1474TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1475TYPE(georef_coord_array) :: this
1476
1477INTEGER :: lsize
1478
1479IF (PRESENT(x) .AND. PRESENT(y)) THEN
1480 lsize = min(SIZE(x), SIZE(y))
1481 ALLOCATE(this%coord(lsize))
1482 this%coord(1:lsize)%x = x(1:lsize)
1483 this%coord(1:lsize)%y = y(1:lsize)
1484ENDIF
1485this%topo = optio_l(topo)
1487
1488END FUNCTION georef_coord_array_new
1489
1490
1491SUBROUTINE georef_coord_array_delete(this)
1492TYPE(georef_coord_array),INTENT(inout) :: this
1493
1494TYPE(georef_coord_array) :: lobj
1495
1496this = lobj
1497
1498END SUBROUTINE georef_coord_array_delete
1499
1500
1501ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1502TYPE(georef_coord_array),INTENT(in) :: this
1503LOGICAL :: res
1504
1505res = ALLOCATED(this%coord)
1506
1507END FUNCTION georef_coord_array_c_e
1508
1509
1514SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1515TYPE(georef_coord_array),INTENT(in) :: this
1516DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1517DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1518! allocatable per vedere di nascosto l'effetto che fa
1519INTEGER,OPTIONAL,INTENT(out) :: topo
1520TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1521
1522
1523IF (PRESENT(x)) THEN
1524 IF (ALLOCATED(this%coord)) THEN
1525 x = this%coord%x
1526 ENDIF
1527ENDIF
1528IF (PRESENT(y)) THEN
1529 IF (ALLOCATED(this%coord)) THEN
1530 y = this%coord%y
1531 ENDIF
1532ENDIF
1533IF (PRESENT(topo)) topo = this%topo
1535
1536END SUBROUTINE georef_coord_array_getval
1537
1538
1544SUBROUTINE georef_coord_array_compute_bbox(this)
1545TYPE(georef_coord_array),INTENT(inout) :: this
1546
1547IF (ALLOCATED(this%coord)) THEN
1548 this%bbox(1)%x = minval(this%coord(:)%x)
1549 this%bbox(1)%y = minval(this%coord(:)%y)
1550 this%bbox(2)%x = maxval(this%coord(:)%x)
1551 this%bbox(2)%y = maxval(this%coord(:)%y)
1552 this%bbox_updated = .true.
1553ENDIF
1554
1555END SUBROUTINE georef_coord_array_compute_bbox
1556
1557#ifdef HAVE_SHAPELIB
1558! internal method for importing a single shape
1559SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1560TYPE(georef_coord_array),INTENT(OUT) :: this
1561TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1562INTEGER,INTENT(IN) :: nshp
1563
1564TYPE(shpobject) :: shpobj
1565
1566! read shape object
1567shpobj = shpreadobject(shphandle, nshp)
1568IF (.NOT.shpisnull(shpobj)) THEN
1569! import it in georef_coord object
1570 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1571 topo=shpobj%nshptype)
1572 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1573 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1574 ELSE IF (ALLOCATED(this%parts)) THEN
1575 DEALLOCATE(this%parts)
1576 ENDIF
1577 CALL shpdestroyobject(shpobj)
1578 CALL compute_bbox(this)
1579ENDIF
1580
1581
1582END SUBROUTINE georef_coord_array_import
1583
1584
1585! internal method for exporting a single shape
1586SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1587TYPE(georef_coord_array),INTENT(in) :: this
1588TYPE(shpfileobject),INTENT(inout) :: shphandle
1589INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1590
1591INTEGER :: i
1592TYPE(shpobject) :: shpobj
1593
1594IF (ALLOCATED(this%coord)) THEN
1595 IF (ALLOCATED(this%parts)) THEN
1596 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1597 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1598 ELSE
1599 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1600 this%coord(:)%x, this%coord(:)%y)
1601 ENDIF
1602ELSE
1603 RETURN
1604ENDIF
1605
1606IF (.NOT.shpisnull(shpobj)) THEN
1607 i = shpwriteobject(shphandle, nshp, shpobj)
1608 CALL shpdestroyobject(shpobj)
1609ENDIF
1610
1611END SUBROUTINE georef_coord_array_export
1612
1613
1624SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1625TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1626CHARACTER(len=*),INTENT(in) :: shpfile
1627
1628REAL(kind=fp_d) :: minb(4), maxb(4)
1629INTEGER :: i, ns, shptype, dbfnf, dbfnr
1630TYPE(shpfileobject) :: shphandle
1631
1632shphandle = shpopen(trim(shpfile), 'rb')
1633IF (shpfileisnull(shphandle)) THEN
1634 ! log here
1635 CALL raise_error()
1636 RETURN
1637ENDIF
1638
1639! get info about file
1640CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1641IF (ns > 0) THEN ! allocate and read the object
1643 DO i = 1, ns
1644 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1645 ENDDO
1646ENDIF
1647
1648CALL shpclose(shphandle)
1649! pack object to save memory
1651
1652END SUBROUTINE arrayof_georef_coord_array_import
1653
1654
1660SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1661TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1662CHARACTER(len=*),INTENT(in) :: shpfile
1663
1664INTEGER :: i
1665TYPE(shpfileobject) :: shphandle
1666
1667IF (this%arraysize > 0) THEN
1668 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1669ELSE
1670 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1671ENDIF
1672IF (shpfileisnull(shphandle)) THEN
1673 ! log here
1674 CALL raise_error()
1675 RETURN
1676ENDIF
1677
1678DO i = 1, this%arraysize
1679 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1680ENDDO
1681
1682CALL shpclose(shphandle)
1683
1684END SUBROUTINE arrayof_georef_coord_array_export
1685#endif
1686
1698FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1699TYPE(georef_coord), INTENT(IN) :: this
1700TYPE(georef_coord_array), INTENT(IN) :: poly
1701LOGICAL :: inside
1702
1703INTEGER :: i
1704
1705inside = .false.
1707IF (.NOT.ALLOCATED(poly%coord)) RETURN
1708! if outside bounding box stop here
1709IF (poly%bbox_updated) THEN
1710 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1711ENDIF
1712
1713IF (ALLOCATED(poly%parts)) THEN
1714 DO i = 1, SIZE(poly%parts)-1
1716 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1717 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1718 ENDDO
1719 IF (SIZE(poly%parts) > 0) THEN ! safety check
1721 poly%coord(poly%parts(i)+1:)%x, &
1722 poly%coord(poly%parts(i)+1:)%y)
1723 ENDIF
1724
1725ELSE
1726 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1727 inside = pointinpoly(this%x, this%y, &
1728 poly%coord(:)%x, poly%coord(:)%y)
1729ENDIF
1730
1731CONTAINS
1732
1733FUNCTION pointinpoly(x, y, px, py)
1734DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1735LOGICAL :: pointinpoly
1736
1737INTEGER :: i, j, starti
1738
1739pointinpoly = .false.
1740
1741IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1742 starti = 2
1743 j = 1
1744ELSE ! unclosed polygon
1745 starti = 1
1746 j = SIZE(px)
1747ENDIF
1748DO i = starti, SIZE(px)
1749 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1750 (py(j) <= y .AND. y < py(i))) THEN
1751 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1752 pointinpoly = .NOT. pointinpoly
1753 ENDIF
1754 ENDIF
1755 j = i
1756ENDDO
1757
1758END FUNCTION pointinpoly
1759
1760END FUNCTION georef_coord_inside
1761
1762
1763
Compute forward coordinate transformation from geographical system to projected system. Definition: geo_proj_class.F90:295 Compute backward coordinate transformation from projected system to geographical system. Definition: geo_proj_class.F90:301 Quick method to append an element to the array. Definition: georef_coord_class.F90:401 Compute the distance in m between two points. Definition: georef_coord_class.F90:344 Export an array of georef_coord_array objects to a file in ESRI/Shapefile format. Definition: georef_coord_class.F90:321 Methods for returning the value of object members. Definition: georef_coord_class.F90:282 Import an array of georef_coord_array objects from a file in ESRI/Shapefile format. Definition: georef_coord_class.F90:315 Method for inserting elements of the array at a desired position. Definition: georef_coord_class.F90:392 Determine whether a point lies inside a polygon or a rectangle. Definition: georef_coord_class.F90:339 Method for packing the array object reducing at a minimum the memory occupation, without destroying i... Definition: georef_coord_class.F90:424 Read a single georef_coord object or an array of georef_coord objects from a Fortran FORMATTED or UNF... Definition: georef_coord_class.F90:328 Method for removing elements of the array at a desired position. Definition: georef_coord_class.F90:407 Write a single georef_coord object or an array of georef_coord objects to a Fortran FORMATTED or UNFO... Definition: georef_coord_class.F90:334 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:227 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:253 Derive type defining a single georeferenced point, either in geodetic or in projected coordinates. Definition: georef_coord_class.F90:241 |