libsim Versione 7.2.1
|
◆ georef_coord_array_getval()
Query a georef_coord_array object. This is the correct way to retrieve the contents of a georef_coord_array object, since its members are declared as PRIVATE.
Definizione alla linea 1005 del file georef_coord_class.F90. 1006! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1007! authors:
1008! Davide Cesari <dcesari@arpa.emr.it>
1009! Paolo Patruno <ppatruno@arpa.emr.it>
1010
1011! This program is free software; you can redistribute it and/or
1012! modify it under the terms of the GNU General Public License as
1013! published by the Free Software Foundation; either version 2 of
1014! the License, or (at your option) any later version.
1015
1016! This program is distributed in the hope that it will be useful,
1017! but WITHOUT ANY WARRANTY; without even the implied warranty of
1018! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1019! GNU General Public License for more details.
1020
1021! You should have received a copy of the GNU General Public License
1022! along with this program. If not, see <http://www.gnu.org/licenses/>.
1023#include "config.h"
1024
1042USE geo_proj_class
1043#ifdef HAVE_SHAPELIB
1044USE shapelib
1045#endif
1046IMPLICIT NONE
1047
1053 PRIVATE
1054 DOUBLE PRECISION :: x=dmiss, y=dmiss
1056
1059
1065 PRIVATE
1066 INTEGER,ALLOCATABLE :: parts(:)
1067 TYPE(georef_coord),ALLOCATABLE :: coord(:)
1068 INTEGER :: topo=imiss
1069 TYPE(geo_proj) :: proj
1070 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
1071 LOGICAL :: bbox_updated=.false.
1073
1074INTEGER,PARAMETER :: georef_coord_array_point = 1
1075INTEGER,PARAMETER :: georef_coord_array_arc = 3
1076INTEGER,PARAMETER :: georef_coord_array_polygon = 5
1077INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
1078
1079
1084 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1085END INTERFACE
1086
1089 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1090END INTERFACE
1091
1094 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
1095END INTERFACE
1096
1097INTERFACE compute_bbox
1098 MODULE PROCEDURE georef_coord_array_compute_bbox
1099END INTERFACE
1100
1102INTERFACE OPERATOR (==)
1103 MODULE PROCEDURE georef_coord_eq
1104END INTERFACE
1105
1107INTERFACE OPERATOR (/=)
1108 MODULE PROCEDURE georef_coord_ne
1109END INTERFACE
1110
1113INTERFACE OPERATOR (>=)
1114 MODULE PROCEDURE georef_coord_ge
1115END INTERFACE
1116
1119INTERFACE OPERATOR (<=)
1120 MODULE PROCEDURE georef_coord_le
1121END INTERFACE
1122
1123#ifdef HAVE_SHAPELIB
1124
1126INTERFACE import
1127 MODULE PROCEDURE arrayof_georef_coord_array_import
1128END INTERFACE
1129
1133 MODULE PROCEDURE arrayof_georef_coord_array_export
1134END INTERFACE
1135#endif
1136
1140 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1141END INTERFACE
1142
1146 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1147END INTERFACE
1148
1151 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1152END INTERFACE
1153
1156 MODULE PROCEDURE georef_coord_dist
1157END INTERFACE
1158
1159#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1160#define ARRAYOF_TYPE arrayof_georef_coord_array
1161!define ARRAYOF_ORIGEQ 0
1162#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1163#include "arrayof_pre.F90"
1164! from arrayof
1166
1167PRIVATE
1169 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1170 georef_coord_array_polygon, georef_coord_array_multipoint, &
1172#ifdef HAVE_SHAPELIB
1174#endif
1176 georef_coord_new, georef_coord_array_new
1177
1178CONTAINS
1179
1180#include "arrayof_post.F90"
1181
1182! ===================
1183! == georef_coord ==
1184! ===================
1188FUNCTION georef_coord_new(x, y) RESULT(this)
1189DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1190DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1191TYPE(georef_coord) :: this
1192
1195
1196END FUNCTION georef_coord_new
1197
1198
1199SUBROUTINE georef_coord_delete(this)
1200TYPE(georef_coord),INTENT(inout) :: this
1201
1202this%x = dmiss
1203this%y = dmiss
1204
1205END SUBROUTINE georef_coord_delete
1206
1207
1208ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1209TYPE(georef_coord),INTENT(in) :: this
1210LOGICAL :: res
1211
1212res = .NOT. this == georef_coord_miss
1213
1214END FUNCTION georef_coord_c_e
1215
1216
1223ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1224TYPE(georef_coord),INTENT(in) :: this
1225DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1226DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1227
1228IF (PRESENT(x)) x = this%x
1229IF (PRESENT(y)) y = this%y
1230
1231END SUBROUTINE georef_coord_getval
1232
1233
1242ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1243TYPE(georef_coord),INTENT(in) :: this
1244TYPE(geo_proj),INTENT(in) :: proj
1245DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1246DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1247DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1248DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1249
1250DOUBLE PRECISION :: llon, llat
1251
1252IF (PRESENT(x)) x = this%x
1253IF (PRESENT(y)) y = this%y
1254IF (PRESENT(lon) .OR. present(lat)) THEN
1256 IF (PRESENT(lon)) lon = llon
1257 IF (PRESENT(lat)) lat = llat
1258ENDIF
1259
1260END SUBROUTINE georef_coord_proj_getval
1261
1262
1263! document and improve
1264ELEMENTAL FUNCTION getlat(this)
1265TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1266DOUBLE PRECISION :: getlat ! latitudine geografica
1267
1268getlat = this%y ! change!!!
1269
1270END FUNCTION getlat
1271
1272! document and improve
1273ELEMENTAL FUNCTION getlon(this)
1274TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1275DOUBLE PRECISION :: getlon ! longitudine geografica
1276
1277getlon = this%x ! change!!!
1278
1279END FUNCTION getlon
1280
1281
1282ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1283TYPE(georef_coord),INTENT(IN) :: this, that
1284LOGICAL :: res
1285
1286res = (this%x == that%x .AND. this%y == that%y)
1287
1288END FUNCTION georef_coord_eq
1289
1290
1291ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1292TYPE(georef_coord),INTENT(IN) :: this, that
1293LOGICAL :: res
1294
1295res = (this%x >= that%x .AND. this%y >= that%y)
1296
1297END FUNCTION georef_coord_ge
1298
1299
1300ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1301TYPE(georef_coord),INTENT(IN) :: this, that
1302LOGICAL :: res
1303
1304res = (this%x <= that%x .AND. this%y <= that%y)
1305
1306END FUNCTION georef_coord_le
1307
1308
1309ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1310TYPE(georef_coord),INTENT(IN) :: this, that
1311LOGICAL :: res
1312
1313res = .NOT.(this == that)
1314
1315END FUNCTION georef_coord_ne
1316
1317
1323SUBROUTINE georef_coord_read_unit(this, unit)
1324TYPE(georef_coord),INTENT(out) :: this
1325INTEGER, INTENT(in) :: unit
1326
1327CALL georef_coord_vect_read_unit((/this/), unit)
1328
1329END SUBROUTINE georef_coord_read_unit
1330
1331
1337SUBROUTINE georef_coord_vect_read_unit(this, unit)
1338TYPE(georef_coord) :: this(:)
1339INTEGER, INTENT(in) :: unit
1340
1341CHARACTER(len=40) :: form
1342INTEGER :: i
1343
1344INQUIRE(unit, form=form)
1345IF (form == 'FORMATTED') THEN
1346 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1347!TODO bug gfortran compiler !
1348!missing values are unredeable when formatted
1349ELSE
1350 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1351ENDIF
1352
1353END SUBROUTINE georef_coord_vect_read_unit
1354
1355
1360SUBROUTINE georef_coord_write_unit(this, unit)
1361TYPE(georef_coord),INTENT(in) :: this
1362INTEGER, INTENT(in) :: unit
1363
1364CALL georef_coord_vect_write_unit((/this/), unit)
1365
1366END SUBROUTINE georef_coord_write_unit
1367
1368
1373SUBROUTINE georef_coord_vect_write_unit(this, unit)
1374TYPE(georef_coord),INTENT(in) :: this(:)
1375INTEGER, INTENT(in) :: unit
1376
1377CHARACTER(len=40) :: form
1378INTEGER :: i
1379
1380INQUIRE(unit, form=form)
1381IF (form == 'FORMATTED') THEN
1382 WRITE(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 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1387ENDIF
1388
1389END SUBROUTINE georef_coord_vect_write_unit
1390
1391
1394FUNCTION georef_coord_dist(this, that) RESULT(dist)
1396TYPE(georef_coord), INTENT (IN) :: this
1397TYPE(georef_coord), INTENT (IN) :: that
1398DOUBLE PRECISION :: dist
1399
1400DOUBLE PRECISION :: x,y
1401! Distanza approssimata, valida per piccoli angoli
1402
1403x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1404y = (this%y-that%y)
1405dist = sqrt(x**2 + y**2)*degrad*rearth
1406
1407END FUNCTION georef_coord_dist
1408
1409
1415FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1416TYPE(georef_coord),INTENT(IN) :: this
1417TYPE(georef_coord),INTENT(IN) :: coordmin
1418TYPE(georef_coord),INTENT(IN) :: coordmax
1419LOGICAL :: res
1420
1421res = (this >= coordmin .AND. this <= coordmax)
1422
1423END FUNCTION georef_coord_inside_rectang
1424
1425
1426! ========================
1427! == georef_coord_array ==
1428! ========================
1434FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1435DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1436DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1437INTEGER,INTENT(in),OPTIONAL :: topo
1438TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1439TYPE(georef_coord_array) :: this
1440
1441INTEGER :: lsize
1442
1443IF (PRESENT(x) .AND. PRESENT(y)) THEN
1444 lsize = min(SIZE(x), SIZE(y))
1445 ALLOCATE(this%coord(lsize))
1446 this%coord(1:lsize)%x = x(1:lsize)
1447 this%coord(1:lsize)%y = y(1:lsize)
1448ENDIF
1449this%topo = optio_l(topo)
1451
1452END FUNCTION georef_coord_array_new
1453
1454
1455SUBROUTINE georef_coord_array_delete(this)
1456TYPE(georef_coord_array),INTENT(inout) :: this
1457
1458TYPE(georef_coord_array) :: lobj
1459
1460this = lobj
1461
1462END SUBROUTINE georef_coord_array_delete
1463
1464
1465ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1466TYPE(georef_coord_array),INTENT(in) :: this
1467LOGICAL :: res
1468
1469res = ALLOCATED(this%coord)
1470
1471END FUNCTION georef_coord_array_c_e
1472
1473
1478SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1479TYPE(georef_coord_array),INTENT(in) :: this
1480DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1481DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1482! allocatable per vedere di nascosto l'effetto che fa
1483INTEGER,OPTIONAL,INTENT(out) :: topo
1484TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1485
1486
1487IF (PRESENT(x)) THEN
1488 IF (ALLOCATED(this%coord)) THEN
1489 x = this%coord%x
1490 ENDIF
1491ENDIF
1492IF (PRESENT(y)) THEN
1493 IF (ALLOCATED(this%coord)) THEN
1494 y = this%coord%y
1495 ENDIF
1496ENDIF
1497IF (PRESENT(topo)) topo = this%topo
1499
1500END SUBROUTINE georef_coord_array_getval
1501
1502
1508SUBROUTINE georef_coord_array_compute_bbox(this)
1509TYPE(georef_coord_array),INTENT(inout) :: this
1510
1511IF (ALLOCATED(this%coord)) THEN
1512 this%bbox(1)%x = minval(this%coord(:)%x)
1513 this%bbox(1)%y = minval(this%coord(:)%y)
1514 this%bbox(2)%x = maxval(this%coord(:)%x)
1515 this%bbox(2)%y = maxval(this%coord(:)%y)
1516 this%bbox_updated = .true.
1517ENDIF
1518
1519END SUBROUTINE georef_coord_array_compute_bbox
1520
1521#ifdef HAVE_SHAPELIB
1522! internal method for importing a single shape
1523SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1524TYPE(georef_coord_array),INTENT(OUT) :: this
1525TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1526INTEGER,INTENT(IN) :: nshp
1527
1528TYPE(shpobject) :: shpobj
1529
1530! read shape object
1531shpobj = shpreadobject(shphandle, nshp)
1532IF (.NOT.shpisnull(shpobj)) THEN
1533! import it in georef_coord object
1534 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1535 topo=shpobj%nshptype)
1536 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1537 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1538 ELSE IF (ALLOCATED(this%parts)) THEN
1539 DEALLOCATE(this%parts)
1540 ENDIF
1541 CALL shpdestroyobject(shpobj)
1542 CALL compute_bbox(this)
1543ENDIF
1544
1545
1546END SUBROUTINE georef_coord_array_import
1547
1548
1549! internal method for exporting a single shape
1550SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1551TYPE(georef_coord_array),INTENT(in) :: this
1552TYPE(shpfileobject),INTENT(inout) :: shphandle
1553INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1554
1555INTEGER :: i
1556TYPE(shpobject) :: shpobj
1557
1558IF (ALLOCATED(this%coord)) THEN
1559 IF (ALLOCATED(this%parts)) THEN
1560 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1561 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1562 ELSE
1563 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1564 this%coord(:)%x, this%coord(:)%y)
1565 ENDIF
1566ELSE
1567 RETURN
1568ENDIF
1569
1570IF (.NOT.shpisnull(shpobj)) THEN
1571 i = shpwriteobject(shphandle, nshp, shpobj)
1572 CALL shpdestroyobject(shpobj)
1573ENDIF
1574
1575END SUBROUTINE georef_coord_array_export
1576
1577
1588SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1589TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1590CHARACTER(len=*),INTENT(in) :: shpfile
1591
1592REAL(kind=fp_d) :: minb(4), maxb(4)
1593INTEGER :: i, ns, shptype, dbfnf, dbfnr
1594TYPE(shpfileobject) :: shphandle
1595
1596shphandle = shpopen(trim(shpfile), 'rb')
1597IF (shpfileisnull(shphandle)) THEN
1598 ! log here
1599 CALL raise_error()
1600 RETURN
1601ENDIF
1602
1603! get info about file
1604CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1605IF (ns > 0) THEN ! allocate and read the object
1607 DO i = 1, ns
1608 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1609 ENDDO
1610ENDIF
1611
1612CALL shpclose(shphandle)
1613! pack object to save memory
1615
1616END SUBROUTINE arrayof_georef_coord_array_import
1617
1618
1624SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1625TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1626CHARACTER(len=*),INTENT(in) :: shpfile
1627
1628INTEGER :: i
1629TYPE(shpfileobject) :: shphandle
1630
1631IF (this%arraysize > 0) THEN
1632 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1633ELSE
1634 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1635ENDIF
1636IF (shpfileisnull(shphandle)) THEN
1637 ! log here
1638 CALL raise_error()
1639 RETURN
1640ENDIF
1641
1642DO i = 1, this%arraysize
1643 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1644ENDDO
1645
1646CALL shpclose(shphandle)
1647
1648END SUBROUTINE arrayof_georef_coord_array_export
1649#endif
1650
1662FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1663TYPE(georef_coord), INTENT(IN) :: this
1664TYPE(georef_coord_array), INTENT(IN) :: poly
1665LOGICAL :: inside
1666
1667INTEGER :: i
1668
1669inside = .false.
1671IF (.NOT.ALLOCATED(poly%coord)) RETURN
1672! if outside bounding box stop here
1673IF (poly%bbox_updated) THEN
1674 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1675ENDIF
1676
1677IF (ALLOCATED(poly%parts)) THEN
1678 DO i = 1, SIZE(poly%parts)-1
1680 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1681 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1682 ENDDO
1683 IF (SIZE(poly%parts) > 0) THEN ! safety check
1685 poly%coord(poly%parts(i)+1:)%x, &
1686 poly%coord(poly%parts(i)+1:)%y)
1687 ENDIF
1688
1689ELSE
1690 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1691 inside = pointinpoly(this%x, this%y, &
1692 poly%coord(:)%x, poly%coord(:)%y)
1693ENDIF
1694
1695CONTAINS
1696
1697FUNCTION pointinpoly(x, y, px, py)
1698DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1699LOGICAL :: pointinpoly
1700
1701INTEGER :: i, j, starti
1702
1703pointinpoly = .false.
1704
1705IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1706 starti = 2
1707 j = 1
1708ELSE ! unclosed polygon
1709 starti = 1
1710 j = SIZE(px)
1711ENDIF
1712DO i = starti, SIZE(px)
1713 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1714 (py(j) <= y .AND. y < py(i))) THEN
1715 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1716 pointinpoly = .NOT. pointinpoly
1717 ENDIF
1718 ENDIF
1719 j = i
1720ENDDO
1721
1722END FUNCTION pointinpoly
1723
1724END FUNCTION georef_coord_inside
1725
1726
1727
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 |