libsim Versione 7.2.4
|
◆ georef_coord_array_new()
Construct a georef_coord_array object with the optional parameters provided. If coordinates are not provided the object obtained is empty (missing, see c_e function), if coordinate arrays are of different lengths the georef_coord_array is initialised to the shortest length provided.
Definizione alla linea 961 del file georef_coord_class.F90. 962! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
963! authors:
964! Davide Cesari <dcesari@arpa.emr.it>
965! Paolo Patruno <ppatruno@arpa.emr.it>
966
967! This program is free software; you can redistribute it and/or
968! modify it under the terms of the GNU General Public License as
969! published by the Free Software Foundation; either version 2 of
970! the License, or (at your option) any later version.
971
972! This program is distributed in the hope that it will be useful,
973! but WITHOUT ANY WARRANTY; without even the implied warranty of
974! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
975! GNU General Public License for more details.
976
977! You should have received a copy of the GNU General Public License
978! along with this program. If not, see <http://www.gnu.org/licenses/>.
979#include "config.h"
980
998USE geo_proj_class
999#ifdef HAVE_SHAPELIB
1000USE shapelib
1001#endif
1002IMPLICIT NONE
1003
1009 PRIVATE
1010 DOUBLE PRECISION :: x=dmiss, y=dmiss
1012
1015
1021 PRIVATE
1022 INTEGER,ALLOCATABLE :: parts(:)
1023 TYPE(georef_coord),ALLOCATABLE :: coord(:)
1024 INTEGER :: topo=imiss
1025 TYPE(geo_proj) :: proj
1026 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
1027 LOGICAL :: bbox_updated=.false.
1029
1030INTEGER,PARAMETER :: georef_coord_array_point = 1
1031INTEGER,PARAMETER :: georef_coord_array_arc = 3
1032INTEGER,PARAMETER :: georef_coord_array_polygon = 5
1033INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
1034
1035
1040 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1041END INTERFACE
1042
1045 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1046END INTERFACE
1047
1050 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
1051END INTERFACE
1052
1053INTERFACE compute_bbox
1054 MODULE PROCEDURE georef_coord_array_compute_bbox
1055END INTERFACE
1056
1058INTERFACE OPERATOR (==)
1059 MODULE PROCEDURE georef_coord_eq
1060END INTERFACE
1061
1063INTERFACE OPERATOR (/=)
1064 MODULE PROCEDURE georef_coord_ne
1065END INTERFACE
1066
1069INTERFACE OPERATOR (>=)
1070 MODULE PROCEDURE georef_coord_ge
1071END INTERFACE
1072
1075INTERFACE OPERATOR (<=)
1076 MODULE PROCEDURE georef_coord_le
1077END INTERFACE
1078
1079#ifdef HAVE_SHAPELIB
1080
1082INTERFACE import
1083 MODULE PROCEDURE arrayof_georef_coord_array_import
1084END INTERFACE
1085
1089 MODULE PROCEDURE arrayof_georef_coord_array_export
1090END INTERFACE
1091#endif
1092
1096 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1097END INTERFACE
1098
1102 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1103END INTERFACE
1104
1107 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1108END INTERFACE
1109
1112 MODULE PROCEDURE georef_coord_dist
1113END INTERFACE
1114
1115#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1116#define ARRAYOF_TYPE arrayof_georef_coord_array
1117!define ARRAYOF_ORIGEQ 0
1118#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1119#include "arrayof_pre.F90"
1120! from arrayof
1122
1123PRIVATE
1125 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1126 georef_coord_array_polygon, georef_coord_array_multipoint, &
1128#ifdef HAVE_SHAPELIB
1130#endif
1132 georef_coord_new, georef_coord_array_new
1133
1134CONTAINS
1135
1136#include "arrayof_post.F90"
1137
1138! ===================
1139! == georef_coord ==
1140! ===================
1144FUNCTION georef_coord_new(x, y) RESULT(this)
1145DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1146DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1147TYPE(georef_coord) :: this
1148
1151
1152END FUNCTION georef_coord_new
1153
1154
1155SUBROUTINE georef_coord_delete(this)
1156TYPE(georef_coord),INTENT(inout) :: this
1157
1158this%x = dmiss
1159this%y = dmiss
1160
1161END SUBROUTINE georef_coord_delete
1162
1163
1164ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1165TYPE(georef_coord),INTENT(in) :: this
1166LOGICAL :: res
1167
1168res = .NOT. this == georef_coord_miss
1169
1170END FUNCTION georef_coord_c_e
1171
1172
1179ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1180TYPE(georef_coord),INTENT(in) :: this
1181DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1182DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1183
1184IF (PRESENT(x)) x = this%x
1185IF (PRESENT(y)) y = this%y
1186
1187END SUBROUTINE georef_coord_getval
1188
1189
1198ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1199TYPE(georef_coord),INTENT(in) :: this
1200TYPE(geo_proj),INTENT(in) :: proj
1201DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1202DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1203DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1204DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1205
1206DOUBLE PRECISION :: llon, llat
1207
1208IF (PRESENT(x)) x = this%x
1209IF (PRESENT(y)) y = this%y
1210IF (PRESENT(lon) .OR. present(lat)) THEN
1212 IF (PRESENT(lon)) lon = llon
1213 IF (PRESENT(lat)) lat = llat
1214ENDIF
1215
1216END SUBROUTINE georef_coord_proj_getval
1217
1218
1219! document and improve
1220ELEMENTAL FUNCTION getlat(this)
1221TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1222DOUBLE PRECISION :: getlat ! latitudine geografica
1223
1224getlat = this%y ! change!!!
1225
1226END FUNCTION getlat
1227
1228! document and improve
1229ELEMENTAL FUNCTION getlon(this)
1230TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1231DOUBLE PRECISION :: getlon ! longitudine geografica
1232
1233getlon = this%x ! change!!!
1234
1235END FUNCTION getlon
1236
1237
1238ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1239TYPE(georef_coord),INTENT(IN) :: this, that
1240LOGICAL :: res
1241
1242res = (this%x == that%x .AND. this%y == that%y)
1243
1244END FUNCTION georef_coord_eq
1245
1246
1247ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1248TYPE(georef_coord),INTENT(IN) :: this, that
1249LOGICAL :: res
1250
1251res = (this%x >= that%x .AND. this%y >= that%y)
1252
1253END FUNCTION georef_coord_ge
1254
1255
1256ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1257TYPE(georef_coord),INTENT(IN) :: this, that
1258LOGICAL :: res
1259
1260res = (this%x <= that%x .AND. this%y <= that%y)
1261
1262END FUNCTION georef_coord_le
1263
1264
1265ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1266TYPE(georef_coord),INTENT(IN) :: this, that
1267LOGICAL :: res
1268
1269res = .NOT.(this == that)
1270
1271END FUNCTION georef_coord_ne
1272
1273
1279SUBROUTINE georef_coord_read_unit(this, unit)
1280TYPE(georef_coord),INTENT(out) :: this
1281INTEGER, INTENT(in) :: unit
1282
1283CALL georef_coord_vect_read_unit((/this/), unit)
1284
1285END SUBROUTINE georef_coord_read_unit
1286
1287
1293SUBROUTINE georef_coord_vect_read_unit(this, unit)
1294TYPE(georef_coord) :: this(:)
1295INTEGER, INTENT(in) :: unit
1296
1297CHARACTER(len=40) :: form
1298INTEGER :: i
1299
1300INQUIRE(unit, form=form)
1301IF (form == 'FORMATTED') THEN
1302 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1303!TODO bug gfortran compiler !
1304!missing values are unredeable when formatted
1305ELSE
1306 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1307ENDIF
1308
1309END SUBROUTINE georef_coord_vect_read_unit
1310
1311
1316SUBROUTINE georef_coord_write_unit(this, unit)
1317TYPE(georef_coord),INTENT(in) :: this
1318INTEGER, INTENT(in) :: unit
1319
1320CALL georef_coord_vect_write_unit((/this/), unit)
1321
1322END SUBROUTINE georef_coord_write_unit
1323
1324
1329SUBROUTINE georef_coord_vect_write_unit(this, unit)
1330TYPE(georef_coord),INTENT(in) :: this(:)
1331INTEGER, INTENT(in) :: unit
1332
1333CHARACTER(len=40) :: form
1334INTEGER :: i
1335
1336INQUIRE(unit, form=form)
1337IF (form == 'FORMATTED') THEN
1338 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1339!TODO bug gfortran compiler !
1340!missing values are unredeable when formatted
1341ELSE
1342 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1343ENDIF
1344
1345END SUBROUTINE georef_coord_vect_write_unit
1346
1347
1350FUNCTION georef_coord_dist(this, that) RESULT(dist)
1352TYPE(georef_coord), INTENT (IN) :: this
1353TYPE(georef_coord), INTENT (IN) :: that
1354DOUBLE PRECISION :: dist
1355
1356DOUBLE PRECISION :: x,y
1357! Distanza approssimata, valida per piccoli angoli
1358
1359x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1360y = (this%y-that%y)
1361dist = sqrt(x**2 + y**2)*degrad*rearth
1362
1363END FUNCTION georef_coord_dist
1364
1365
1371FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1372TYPE(georef_coord),INTENT(IN) :: this
1373TYPE(georef_coord),INTENT(IN) :: coordmin
1374TYPE(georef_coord),INTENT(IN) :: coordmax
1375LOGICAL :: res
1376
1377res = (this >= coordmin .AND. this <= coordmax)
1378
1379END FUNCTION georef_coord_inside_rectang
1380
1381
1382! ========================
1383! == georef_coord_array ==
1384! ========================
1390FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1391DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1392DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1393INTEGER,INTENT(in),OPTIONAL :: topo
1394TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1395TYPE(georef_coord_array) :: this
1396
1397INTEGER :: lsize
1398
1399IF (PRESENT(x) .AND. PRESENT(y)) THEN
1400 lsize = min(SIZE(x), SIZE(y))
1401 ALLOCATE(this%coord(lsize))
1402 this%coord(1:lsize)%x = x(1:lsize)
1403 this%coord(1:lsize)%y = y(1:lsize)
1404ENDIF
1405this%topo = optio_l(topo)
1407
1408END FUNCTION georef_coord_array_new
1409
1410
1411SUBROUTINE georef_coord_array_delete(this)
1412TYPE(georef_coord_array),INTENT(inout) :: this
1413
1414TYPE(georef_coord_array) :: lobj
1415
1416this = lobj
1417
1418END SUBROUTINE georef_coord_array_delete
1419
1420
1421ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1422TYPE(georef_coord_array),INTENT(in) :: this
1423LOGICAL :: res
1424
1425res = ALLOCATED(this%coord)
1426
1427END FUNCTION georef_coord_array_c_e
1428
1429
1434SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1435TYPE(georef_coord_array),INTENT(in) :: this
1436DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1437DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1438! allocatable per vedere di nascosto l'effetto che fa
1439INTEGER,OPTIONAL,INTENT(out) :: topo
1440TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1441
1442
1443IF (PRESENT(x)) THEN
1444 IF (ALLOCATED(this%coord)) THEN
1445 x = this%coord%x
1446 ENDIF
1447ENDIF
1448IF (PRESENT(y)) THEN
1449 IF (ALLOCATED(this%coord)) THEN
1450 y = this%coord%y
1451 ENDIF
1452ENDIF
1453IF (PRESENT(topo)) topo = this%topo
1455
1456END SUBROUTINE georef_coord_array_getval
1457
1458
1464SUBROUTINE georef_coord_array_compute_bbox(this)
1465TYPE(georef_coord_array),INTENT(inout) :: this
1466
1467IF (ALLOCATED(this%coord)) THEN
1468 this%bbox(1)%x = minval(this%coord(:)%x)
1469 this%bbox(1)%y = minval(this%coord(:)%y)
1470 this%bbox(2)%x = maxval(this%coord(:)%x)
1471 this%bbox(2)%y = maxval(this%coord(:)%y)
1472 this%bbox_updated = .true.
1473ENDIF
1474
1475END SUBROUTINE georef_coord_array_compute_bbox
1476
1477#ifdef HAVE_SHAPELIB
1478! internal method for importing a single shape
1479SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1480TYPE(georef_coord_array),INTENT(OUT) :: this
1481TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1482INTEGER,INTENT(IN) :: nshp
1483
1484TYPE(shpobject) :: shpobj
1485
1486! read shape object
1487shpobj = shpreadobject(shphandle, nshp)
1488IF (.NOT.shpisnull(shpobj)) THEN
1489! import it in georef_coord object
1490 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1491 topo=shpobj%nshptype)
1492 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1493 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1494 ELSE IF (ALLOCATED(this%parts)) THEN
1495 DEALLOCATE(this%parts)
1496 ENDIF
1497 CALL shpdestroyobject(shpobj)
1498 CALL compute_bbox(this)
1499ENDIF
1500
1501
1502END SUBROUTINE georef_coord_array_import
1503
1504
1505! internal method for exporting a single shape
1506SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1507TYPE(georef_coord_array),INTENT(in) :: this
1508TYPE(shpfileobject),INTENT(inout) :: shphandle
1509INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1510
1511INTEGER :: i
1512TYPE(shpobject) :: shpobj
1513
1514IF (ALLOCATED(this%coord)) THEN
1515 IF (ALLOCATED(this%parts)) THEN
1516 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1517 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1518 ELSE
1519 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1520 this%coord(:)%x, this%coord(:)%y)
1521 ENDIF
1522ELSE
1523 RETURN
1524ENDIF
1525
1526IF (.NOT.shpisnull(shpobj)) THEN
1527 i = shpwriteobject(shphandle, nshp, shpobj)
1528 CALL shpdestroyobject(shpobj)
1529ENDIF
1530
1531END SUBROUTINE georef_coord_array_export
1532
1533
1544SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1545TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1546CHARACTER(len=*),INTENT(in) :: shpfile
1547
1548REAL(kind=fp_d) :: minb(4), maxb(4)
1549INTEGER :: i, ns, shptype, dbfnf, dbfnr
1550TYPE(shpfileobject) :: shphandle
1551
1552shphandle = shpopen(trim(shpfile), 'rb')
1553IF (shpfileisnull(shphandle)) THEN
1554 ! log here
1555 CALL raise_error()
1556 RETURN
1557ENDIF
1558
1559! get info about file
1560CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1561IF (ns > 0) THEN ! allocate and read the object
1563 DO i = 1, ns
1564 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1565 ENDDO
1566ENDIF
1567
1568CALL shpclose(shphandle)
1569! pack object to save memory
1571
1572END SUBROUTINE arrayof_georef_coord_array_import
1573
1574
1580SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1581TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1582CHARACTER(len=*),INTENT(in) :: shpfile
1583
1584INTEGER :: i
1585TYPE(shpfileobject) :: shphandle
1586
1587IF (this%arraysize > 0) THEN
1588 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1589ELSE
1590 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1591ENDIF
1592IF (shpfileisnull(shphandle)) THEN
1593 ! log here
1594 CALL raise_error()
1595 RETURN
1596ENDIF
1597
1598DO i = 1, this%arraysize
1599 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1600ENDDO
1601
1602CALL shpclose(shphandle)
1603
1604END SUBROUTINE arrayof_georef_coord_array_export
1605#endif
1606
1618FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1619TYPE(georef_coord), INTENT(IN) :: this
1620TYPE(georef_coord_array), INTENT(IN) :: poly
1621LOGICAL :: inside
1622
1623INTEGER :: i
1624
1625inside = .false.
1627IF (.NOT.ALLOCATED(poly%coord)) RETURN
1628! if outside bounding box stop here
1629IF (poly%bbox_updated) THEN
1630 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1631ENDIF
1632
1633IF (ALLOCATED(poly%parts)) THEN
1634 DO i = 1, SIZE(poly%parts)-1
1636 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1637 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1638 ENDDO
1639 IF (SIZE(poly%parts) > 0) THEN ! safety check
1641 poly%coord(poly%parts(i)+1:)%x, &
1642 poly%coord(poly%parts(i)+1:)%y)
1643 ENDIF
1644
1645ELSE
1646 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1647 inside = pointinpoly(this%x, this%y, &
1648 poly%coord(:)%x, poly%coord(:)%y)
1649ENDIF
1650
1651CONTAINS
1652
1653FUNCTION pointinpoly(x, y, px, py)
1654DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1655LOGICAL :: pointinpoly
1656
1657INTEGER :: i, j, starti
1658
1659pointinpoly = .false.
1660
1661IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1662 starti = 2
1663 j = 1
1664ELSE ! unclosed polygon
1665 starti = 1
1666 j = SIZE(px)
1667ENDIF
1668DO i = starti, SIZE(px)
1669 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1670 (py(j) <= y .AND. y < py(i))) THEN
1671 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1672 pointinpoly = .NOT. pointinpoly
1673 ENDIF
1674 ENDIF
1675 j = i
1676ENDDO
1677
1678END FUNCTION pointinpoly
1679
1680END FUNCTION georef_coord_inside
1681
1682
1683
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 |