libsim Versione 7.1.11
|
◆ 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 1011 del file georef_coord_class.F90. 1012! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1013! authors:
1014! Davide Cesari <dcesari@arpa.emr.it>
1015! Paolo Patruno <ppatruno@arpa.emr.it>
1016
1017! This program is free software; you can redistribute it and/or
1018! modify it under the terms of the GNU General Public License as
1019! published by the Free Software Foundation; either version 2 of
1020! the License, or (at your option) any later version.
1021
1022! This program is distributed in the hope that it will be useful,
1023! but WITHOUT ANY WARRANTY; without even the implied warranty of
1024! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1025! GNU General Public License for more details.
1026
1027! You should have received a copy of the GNU General Public License
1028! along with this program. If not, see <http://www.gnu.org/licenses/>.
1029#include "config.h"
1030
1048USE geo_proj_class
1049#ifdef HAVE_SHAPELIB
1050USE shapelib
1051#endif
1052IMPLICIT NONE
1053
1059 PRIVATE
1060 DOUBLE PRECISION :: x=dmiss, y=dmiss
1062
1065
1071 PRIVATE
1072 INTEGER,ALLOCATABLE :: parts(:)
1073 TYPE(georef_coord),ALLOCATABLE :: coord(:)
1074 INTEGER :: topo=imiss
1075 TYPE(geo_proj) :: proj
1076 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
1077 LOGICAL :: bbox_updated=.false.
1079
1080INTEGER,PARAMETER :: georef_coord_array_point = 1
1081INTEGER,PARAMETER :: georef_coord_array_arc = 3
1082INTEGER,PARAMETER :: georef_coord_array_polygon = 5
1083INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
1084
1085
1090 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1091END INTERFACE
1092
1095 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1096END INTERFACE
1097
1100 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
1101END INTERFACE
1102
1103INTERFACE compute_bbox
1104 MODULE PROCEDURE georef_coord_array_compute_bbox
1105END INTERFACE
1106
1108INTERFACE OPERATOR (==)
1109 MODULE PROCEDURE georef_coord_eq
1110END INTERFACE
1111
1113INTERFACE OPERATOR (/=)
1114 MODULE PROCEDURE georef_coord_ne
1115END INTERFACE
1116
1119INTERFACE OPERATOR (>=)
1120 MODULE PROCEDURE georef_coord_ge
1121END INTERFACE
1122
1125INTERFACE OPERATOR (<=)
1126 MODULE PROCEDURE georef_coord_le
1127END INTERFACE
1128
1129#ifdef HAVE_SHAPELIB
1130
1132INTERFACE import
1133 MODULE PROCEDURE arrayof_georef_coord_array_import
1134END INTERFACE
1135
1139 MODULE PROCEDURE arrayof_georef_coord_array_export
1140END INTERFACE
1141#endif
1142
1146 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1147END INTERFACE
1148
1152 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1153END INTERFACE
1154
1157 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1158END INTERFACE
1159
1162 MODULE PROCEDURE georef_coord_dist
1163END INTERFACE
1164
1165#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1166#define ARRAYOF_TYPE arrayof_georef_coord_array
1167!define ARRAYOF_ORIGEQ 0
1168#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1169#include "arrayof_pre.F90"
1170! from arrayof
1172
1173PRIVATE
1175 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1176 georef_coord_array_polygon, georef_coord_array_multipoint, &
1178#ifdef HAVE_SHAPELIB
1180#endif
1182 georef_coord_new, georef_coord_array_new
1183
1184CONTAINS
1185
1186#include "arrayof_post.F90"
1187
1188! ===================
1189! == georef_coord ==
1190! ===================
1194FUNCTION georef_coord_new(x, y) RESULT(this)
1195DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1196DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1197TYPE(georef_coord) :: this
1198
1201
1202END FUNCTION georef_coord_new
1203
1204
1205SUBROUTINE georef_coord_delete(this)
1206TYPE(georef_coord),INTENT(inout) :: this
1207
1208this%x = dmiss
1209this%y = dmiss
1210
1211END SUBROUTINE georef_coord_delete
1212
1213
1214ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1215TYPE(georef_coord),INTENT(in) :: this
1216LOGICAL :: res
1217
1218res = .NOT. this == georef_coord_miss
1219
1220END FUNCTION georef_coord_c_e
1221
1222
1229ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1230TYPE(georef_coord),INTENT(in) :: this
1231DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1232DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1233
1234IF (PRESENT(x)) x = this%x
1235IF (PRESENT(y)) y = this%y
1236
1237END SUBROUTINE georef_coord_getval
1238
1239
1248ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1249TYPE(georef_coord),INTENT(in) :: this
1250TYPE(geo_proj),INTENT(in) :: proj
1251DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1252DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1253DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1254DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1255
1256DOUBLE PRECISION :: llon, llat
1257
1258IF (PRESENT(x)) x = this%x
1259IF (PRESENT(y)) y = this%y
1260IF (PRESENT(lon) .OR. present(lat)) THEN
1262 IF (PRESENT(lon)) lon = llon
1263 IF (PRESENT(lat)) lat = llat
1264ENDIF
1265
1266END SUBROUTINE georef_coord_proj_getval
1267
1268
1269! document and improve
1270ELEMENTAL FUNCTION getlat(this)
1271TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1272DOUBLE PRECISION :: getlat ! latitudine geografica
1273
1274getlat = this%y ! change!!!
1275
1276END FUNCTION getlat
1277
1278! document and improve
1279ELEMENTAL FUNCTION getlon(this)
1280TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1281DOUBLE PRECISION :: getlon ! longitudine geografica
1282
1283getlon = this%x ! change!!!
1284
1285END FUNCTION getlon
1286
1287
1288ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1289TYPE(georef_coord),INTENT(IN) :: this, that
1290LOGICAL :: res
1291
1292res = (this%x == that%x .AND. this%y == that%y)
1293
1294END FUNCTION georef_coord_eq
1295
1296
1297ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1298TYPE(georef_coord),INTENT(IN) :: this, that
1299LOGICAL :: res
1300
1301res = (this%x >= that%x .AND. this%y >= that%y)
1302
1303END FUNCTION georef_coord_ge
1304
1305
1306ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1307TYPE(georef_coord),INTENT(IN) :: this, that
1308LOGICAL :: res
1309
1310res = (this%x <= that%x .AND. this%y <= that%y)
1311
1312END FUNCTION georef_coord_le
1313
1314
1315ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1316TYPE(georef_coord),INTENT(IN) :: this, that
1317LOGICAL :: res
1318
1319res = .NOT.(this == that)
1320
1321END FUNCTION georef_coord_ne
1322
1323
1329SUBROUTINE georef_coord_read_unit(this, unit)
1330TYPE(georef_coord),INTENT(out) :: this
1331INTEGER, INTENT(in) :: unit
1332
1333CALL georef_coord_vect_read_unit((/this/), unit)
1334
1335END SUBROUTINE georef_coord_read_unit
1336
1337
1343SUBROUTINE georef_coord_vect_read_unit(this, unit)
1344TYPE(georef_coord) :: this(:)
1345INTEGER, INTENT(in) :: unit
1346
1347CHARACTER(len=40) :: form
1348INTEGER :: i
1349
1350INQUIRE(unit, form=form)
1351IF (form == 'FORMATTED') THEN
1352 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1353!TODO bug gfortran compiler !
1354!missing values are unredeable when formatted
1355ELSE
1356 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1357ENDIF
1358
1359END SUBROUTINE georef_coord_vect_read_unit
1360
1361
1366SUBROUTINE georef_coord_write_unit(this, unit)
1367TYPE(georef_coord),INTENT(in) :: this
1368INTEGER, INTENT(in) :: unit
1369
1370CALL georef_coord_vect_write_unit((/this/), unit)
1371
1372END SUBROUTINE georef_coord_write_unit
1373
1374
1379SUBROUTINE georef_coord_vect_write_unit(this, unit)
1380TYPE(georef_coord),INTENT(in) :: this(:)
1381INTEGER, INTENT(in) :: unit
1382
1383CHARACTER(len=40) :: form
1384INTEGER :: i
1385
1386INQUIRE(unit, form=form)
1387IF (form == 'FORMATTED') THEN
1388 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1389!TODO bug gfortran compiler !
1390!missing values are unredeable when formatted
1391ELSE
1392 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1393ENDIF
1394
1395END SUBROUTINE georef_coord_vect_write_unit
1396
1397
1400FUNCTION georef_coord_dist(this, that) RESULT(dist)
1402TYPE(georef_coord), INTENT (IN) :: this
1403TYPE(georef_coord), INTENT (IN) :: that
1404DOUBLE PRECISION :: dist
1405
1406DOUBLE PRECISION :: x,y
1407! Distanza approssimata, valida per piccoli angoli
1408
1409x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1410y = (this%y-that%y)
1411dist = sqrt(x**2 + y**2)*degrad*rearth
1412
1413END FUNCTION georef_coord_dist
1414
1415
1421FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1422TYPE(georef_coord),INTENT(IN) :: this
1423TYPE(georef_coord),INTENT(IN) :: coordmin
1424TYPE(georef_coord),INTENT(IN) :: coordmax
1425LOGICAL :: res
1426
1427res = (this >= coordmin .AND. this <= coordmax)
1428
1429END FUNCTION georef_coord_inside_rectang
1430
1431
1432! ========================
1433! == georef_coord_array ==
1434! ========================
1440FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1441DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1442DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1443INTEGER,INTENT(in),OPTIONAL :: topo
1444TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1445TYPE(georef_coord_array) :: this
1446
1447INTEGER :: lsize
1448
1449IF (PRESENT(x) .AND. PRESENT(y)) THEN
1450 lsize = min(SIZE(x), SIZE(y))
1451 ALLOCATE(this%coord(lsize))
1452 this%coord(1:lsize)%x = x(1:lsize)
1453 this%coord(1:lsize)%y = y(1:lsize)
1454ENDIF
1455this%topo = optio_l(topo)
1457
1458END FUNCTION georef_coord_array_new
1459
1460
1461SUBROUTINE georef_coord_array_delete(this)
1462TYPE(georef_coord_array),INTENT(inout) :: this
1463
1464TYPE(georef_coord_array) :: lobj
1465
1466this = lobj
1467
1468END SUBROUTINE georef_coord_array_delete
1469
1470
1471ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1472TYPE(georef_coord_array),INTENT(in) :: this
1473LOGICAL :: res
1474
1475res = ALLOCATED(this%coord)
1476
1477END FUNCTION georef_coord_array_c_e
1478
1479
1484SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1485TYPE(georef_coord_array),INTENT(in) :: this
1486DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1487DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1488! allocatable per vedere di nascosto l'effetto che fa
1489INTEGER,OPTIONAL,INTENT(out) :: topo
1490TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1491
1492
1493IF (PRESENT(x)) THEN
1494 IF (ALLOCATED(this%coord)) THEN
1495 x = this%coord%x
1496 ENDIF
1497ENDIF
1498IF (PRESENT(y)) THEN
1499 IF (ALLOCATED(this%coord)) THEN
1500 y = this%coord%y
1501 ENDIF
1502ENDIF
1503IF (PRESENT(topo)) topo = this%topo
1505
1506END SUBROUTINE georef_coord_array_getval
1507
1508
1514SUBROUTINE georef_coord_array_compute_bbox(this)
1515TYPE(georef_coord_array),INTENT(inout) :: this
1516
1517IF (ALLOCATED(this%coord)) THEN
1518 this%bbox(1)%x = minval(this%coord(:)%x)
1519 this%bbox(1)%y = minval(this%coord(:)%y)
1520 this%bbox(2)%x = maxval(this%coord(:)%x)
1521 this%bbox(2)%y = maxval(this%coord(:)%y)
1522 this%bbox_updated = .true.
1523ENDIF
1524
1525END SUBROUTINE georef_coord_array_compute_bbox
1526
1527#ifdef HAVE_SHAPELIB
1528! internal method for importing a single shape
1529SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1530TYPE(georef_coord_array),INTENT(OUT) :: this
1531TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1532INTEGER,INTENT(IN) :: nshp
1533
1534TYPE(shpobject) :: shpobj
1535
1536! read shape object
1537shpobj = shpreadobject(shphandle, nshp)
1538IF (.NOT.shpisnull(shpobj)) THEN
1539! import it in georef_coord object
1540 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1541 topo=shpobj%nshptype)
1542 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1543 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1544 ELSE IF (ALLOCATED(this%parts)) THEN
1545 DEALLOCATE(this%parts)
1546 ENDIF
1547 CALL shpdestroyobject(shpobj)
1548 CALL compute_bbox(this)
1549ENDIF
1550
1551
1552END SUBROUTINE georef_coord_array_import
1553
1554
1555! internal method for exporting a single shape
1556SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1557TYPE(georef_coord_array),INTENT(in) :: this
1558TYPE(shpfileobject),INTENT(inout) :: shphandle
1559INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1560
1561INTEGER :: i
1562TYPE(shpobject) :: shpobj
1563
1564IF (ALLOCATED(this%coord)) THEN
1565 IF (ALLOCATED(this%parts)) THEN
1566 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1567 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1568 ELSE
1569 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1570 this%coord(:)%x, this%coord(:)%y)
1571 ENDIF
1572ELSE
1573 RETURN
1574ENDIF
1575
1576IF (.NOT.shpisnull(shpobj)) THEN
1577 i = shpwriteobject(shphandle, nshp, shpobj)
1578 CALL shpdestroyobject(shpobj)
1579ENDIF
1580
1581END SUBROUTINE georef_coord_array_export
1582
1583
1594SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1595TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1596CHARACTER(len=*),INTENT(in) :: shpfile
1597
1598REAL(kind=fp_d) :: minb(4), maxb(4)
1599INTEGER :: i, ns, shptype, dbfnf, dbfnr
1600TYPE(shpfileobject) :: shphandle
1601
1602shphandle = shpopen(trim(shpfile), 'rb')
1603IF (shpfileisnull(shphandle)) THEN
1604 ! log here
1605 CALL raise_error()
1606 RETURN
1607ENDIF
1608
1609! get info about file
1610CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1611IF (ns > 0) THEN ! allocate and read the object
1613 DO i = 1, ns
1614 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1615 ENDDO
1616ENDIF
1617
1618CALL shpclose(shphandle)
1619! pack object to save memory
1621
1622END SUBROUTINE arrayof_georef_coord_array_import
1623
1624
1630SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1631TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1632CHARACTER(len=*),INTENT(in) :: shpfile
1633
1634INTEGER :: i
1635TYPE(shpfileobject) :: shphandle
1636
1637IF (this%arraysize > 0) THEN
1638 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1639ELSE
1640 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1641ENDIF
1642IF (shpfileisnull(shphandle)) THEN
1643 ! log here
1644 CALL raise_error()
1645 RETURN
1646ENDIF
1647
1648DO i = 1, this%arraysize
1649 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1650ENDDO
1651
1652CALL shpclose(shphandle)
1653
1654END SUBROUTINE arrayof_georef_coord_array_export
1655#endif
1656
1668FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1669TYPE(georef_coord), INTENT(IN) :: this
1670TYPE(georef_coord_array), INTENT(IN) :: poly
1671LOGICAL :: inside
1672
1673INTEGER :: i
1674
1675inside = .false.
1677IF (.NOT.ALLOCATED(poly%coord)) RETURN
1678! if outside bounding box stop here
1679IF (poly%bbox_updated) THEN
1680 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1681ENDIF
1682
1683IF (ALLOCATED(poly%parts)) THEN
1684 DO i = 1, SIZE(poly%parts)-1
1686 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1687 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1688 ENDDO
1689 IF (SIZE(poly%parts) > 0) THEN ! safety check
1691 poly%coord(poly%parts(i)+1:)%x, &
1692 poly%coord(poly%parts(i)+1:)%y)
1693 ENDIF
1694
1695ELSE
1696 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1697 inside = pointinpoly(this%x, this%y, &
1698 poly%coord(:)%x, poly%coord(:)%y)
1699ENDIF
1700
1701CONTAINS
1702
1703FUNCTION pointinpoly(x, y, px, py)
1704DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1705LOGICAL :: pointinpoly
1706
1707INTEGER :: i, j, starti
1708
1709pointinpoly = .false.
1710
1711IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1712 starti = 2
1713 j = 1
1714ELSE ! unclosed polygon
1715 starti = 1
1716 j = SIZE(px)
1717ENDIF
1718DO i = starti, SIZE(px)
1719 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1720 (py(j) <= y .AND. y < py(i))) THEN
1721 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1722 pointinpoly = .NOT. pointinpoly
1723 ENDIF
1724 ENDIF
1725 j = i
1726ENDDO
1727
1728END FUNCTION pointinpoly
1729
1730END FUNCTION georef_coord_inside
1731
1732
1733
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 |