libsim Versione 7.2.1
|
◆ georef_coord_dist()
Restituisce la distanza in m tra 2 oggetti georef_coord. La distanza è calcolata approssimativamente ed è valida per piccoli angoli.
Definizione alla linea 921 del file georef_coord_class.F90. 922! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
923! authors:
924! Davide Cesari <dcesari@arpa.emr.it>
925! Paolo Patruno <ppatruno@arpa.emr.it>
926
927! This program is free software; you can redistribute it and/or
928! modify it under the terms of the GNU General Public License as
929! published by the Free Software Foundation; either version 2 of
930! the License, or (at your option) any later version.
931
932! This program is distributed in the hope that it will be useful,
933! but WITHOUT ANY WARRANTY; without even the implied warranty of
934! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
935! GNU General Public License for more details.
936
937! You should have received a copy of the GNU General Public License
938! along with this program. If not, see <http://www.gnu.org/licenses/>.
939#include "config.h"
940
958USE geo_proj_class
959#ifdef HAVE_SHAPELIB
960USE shapelib
961#endif
962IMPLICIT NONE
963
969 PRIVATE
970 DOUBLE PRECISION :: x=dmiss, y=dmiss
972
975
981 PRIVATE
982 INTEGER,ALLOCATABLE :: parts(:)
983 TYPE(georef_coord),ALLOCATABLE :: coord(:)
984 INTEGER :: topo=imiss
985 TYPE(geo_proj) :: proj
986 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
987 LOGICAL :: bbox_updated=.false.
989
990INTEGER,PARAMETER :: georef_coord_array_point = 1
991INTEGER,PARAMETER :: georef_coord_array_arc = 3
992INTEGER,PARAMETER :: georef_coord_array_polygon = 5
993INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
994
995
1000 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1001END INTERFACE
1002
1005 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1006END INTERFACE
1007
1010 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
1011END INTERFACE
1012
1013INTERFACE compute_bbox
1014 MODULE PROCEDURE georef_coord_array_compute_bbox
1015END INTERFACE
1016
1018INTERFACE OPERATOR (==)
1019 MODULE PROCEDURE georef_coord_eq
1020END INTERFACE
1021
1023INTERFACE OPERATOR (/=)
1024 MODULE PROCEDURE georef_coord_ne
1025END INTERFACE
1026
1029INTERFACE OPERATOR (>=)
1030 MODULE PROCEDURE georef_coord_ge
1031END INTERFACE
1032
1035INTERFACE OPERATOR (<=)
1036 MODULE PROCEDURE georef_coord_le
1037END INTERFACE
1038
1039#ifdef HAVE_SHAPELIB
1040
1042INTERFACE import
1043 MODULE PROCEDURE arrayof_georef_coord_array_import
1044END INTERFACE
1045
1049 MODULE PROCEDURE arrayof_georef_coord_array_export
1050END INTERFACE
1051#endif
1052
1056 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1057END INTERFACE
1058
1062 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1063END INTERFACE
1064
1067 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1068END INTERFACE
1069
1072 MODULE PROCEDURE georef_coord_dist
1073END INTERFACE
1074
1075#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1076#define ARRAYOF_TYPE arrayof_georef_coord_array
1077!define ARRAYOF_ORIGEQ 0
1078#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1079#include "arrayof_pre.F90"
1080! from arrayof
1082
1083PRIVATE
1085 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1086 georef_coord_array_polygon, georef_coord_array_multipoint, &
1088#ifdef HAVE_SHAPELIB
1090#endif
1092 georef_coord_new, georef_coord_array_new
1093
1094CONTAINS
1095
1096#include "arrayof_post.F90"
1097
1098! ===================
1099! == georef_coord ==
1100! ===================
1104FUNCTION georef_coord_new(x, y) RESULT(this)
1105DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1106DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1107TYPE(georef_coord) :: this
1108
1111
1112END FUNCTION georef_coord_new
1113
1114
1115SUBROUTINE georef_coord_delete(this)
1116TYPE(georef_coord),INTENT(inout) :: this
1117
1118this%x = dmiss
1119this%y = dmiss
1120
1121END SUBROUTINE georef_coord_delete
1122
1123
1124ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1125TYPE(georef_coord),INTENT(in) :: this
1126LOGICAL :: res
1127
1128res = .NOT. this == georef_coord_miss
1129
1130END FUNCTION georef_coord_c_e
1131
1132
1139ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1140TYPE(georef_coord),INTENT(in) :: this
1141DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1142DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1143
1144IF (PRESENT(x)) x = this%x
1145IF (PRESENT(y)) y = this%y
1146
1147END SUBROUTINE georef_coord_getval
1148
1149
1158ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1159TYPE(georef_coord),INTENT(in) :: this
1160TYPE(geo_proj),INTENT(in) :: proj
1161DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1162DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1163DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1164DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1165
1166DOUBLE PRECISION :: llon, llat
1167
1168IF (PRESENT(x)) x = this%x
1169IF (PRESENT(y)) y = this%y
1170IF (PRESENT(lon) .OR. present(lat)) THEN
1172 IF (PRESENT(lon)) lon = llon
1173 IF (PRESENT(lat)) lat = llat
1174ENDIF
1175
1176END SUBROUTINE georef_coord_proj_getval
1177
1178
1179! document and improve
1180ELEMENTAL FUNCTION getlat(this)
1181TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1182DOUBLE PRECISION :: getlat ! latitudine geografica
1183
1184getlat = this%y ! change!!!
1185
1186END FUNCTION getlat
1187
1188! document and improve
1189ELEMENTAL FUNCTION getlon(this)
1190TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1191DOUBLE PRECISION :: getlon ! longitudine geografica
1192
1193getlon = this%x ! change!!!
1194
1195END FUNCTION getlon
1196
1197
1198ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1199TYPE(georef_coord),INTENT(IN) :: this, that
1200LOGICAL :: res
1201
1202res = (this%x == that%x .AND. this%y == that%y)
1203
1204END FUNCTION georef_coord_eq
1205
1206
1207ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1208TYPE(georef_coord),INTENT(IN) :: this, that
1209LOGICAL :: res
1210
1211res = (this%x >= that%x .AND. this%y >= that%y)
1212
1213END FUNCTION georef_coord_ge
1214
1215
1216ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1217TYPE(georef_coord),INTENT(IN) :: this, that
1218LOGICAL :: res
1219
1220res = (this%x <= that%x .AND. this%y <= that%y)
1221
1222END FUNCTION georef_coord_le
1223
1224
1225ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1226TYPE(georef_coord),INTENT(IN) :: this, that
1227LOGICAL :: res
1228
1229res = .NOT.(this == that)
1230
1231END FUNCTION georef_coord_ne
1232
1233
1239SUBROUTINE georef_coord_read_unit(this, unit)
1240TYPE(georef_coord),INTENT(out) :: this
1241INTEGER, INTENT(in) :: unit
1242
1243CALL georef_coord_vect_read_unit((/this/), unit)
1244
1245END SUBROUTINE georef_coord_read_unit
1246
1247
1253SUBROUTINE georef_coord_vect_read_unit(this, unit)
1254TYPE(georef_coord) :: this(:)
1255INTEGER, INTENT(in) :: unit
1256
1257CHARACTER(len=40) :: form
1258INTEGER :: i
1259
1260INQUIRE(unit, form=form)
1261IF (form == 'FORMATTED') THEN
1262 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1263!TODO bug gfortran compiler !
1264!missing values are unredeable when formatted
1265ELSE
1266 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1267ENDIF
1268
1269END SUBROUTINE georef_coord_vect_read_unit
1270
1271
1276SUBROUTINE georef_coord_write_unit(this, unit)
1277TYPE(georef_coord),INTENT(in) :: this
1278INTEGER, INTENT(in) :: unit
1279
1280CALL georef_coord_vect_write_unit((/this/), unit)
1281
1282END SUBROUTINE georef_coord_write_unit
1283
1284
1289SUBROUTINE georef_coord_vect_write_unit(this, unit)
1290TYPE(georef_coord),INTENT(in) :: this(:)
1291INTEGER, INTENT(in) :: unit
1292
1293CHARACTER(len=40) :: form
1294INTEGER :: i
1295
1296INQUIRE(unit, form=form)
1297IF (form == 'FORMATTED') THEN
1298 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1299!TODO bug gfortran compiler !
1300!missing values are unredeable when formatted
1301ELSE
1302 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1303ENDIF
1304
1305END SUBROUTINE georef_coord_vect_write_unit
1306
1307
1310FUNCTION georef_coord_dist(this, that) RESULT(dist)
1312TYPE(georef_coord), INTENT (IN) :: this
1313TYPE(georef_coord), INTENT (IN) :: that
1314DOUBLE PRECISION :: dist
1315
1316DOUBLE PRECISION :: x,y
1317! Distanza approssimata, valida per piccoli angoli
1318
1319x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1320y = (this%y-that%y)
1321dist = sqrt(x**2 + y**2)*degrad*rearth
1322
1323END FUNCTION georef_coord_dist
1324
1325
1331FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1332TYPE(georef_coord),INTENT(IN) :: this
1333TYPE(georef_coord),INTENT(IN) :: coordmin
1334TYPE(georef_coord),INTENT(IN) :: coordmax
1335LOGICAL :: res
1336
1337res = (this >= coordmin .AND. this <= coordmax)
1338
1339END FUNCTION georef_coord_inside_rectang
1340
1341
1342! ========================
1343! == georef_coord_array ==
1344! ========================
1350FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1351DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1352DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1353INTEGER,INTENT(in),OPTIONAL :: topo
1354TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1355TYPE(georef_coord_array) :: this
1356
1357INTEGER :: lsize
1358
1359IF (PRESENT(x) .AND. PRESENT(y)) THEN
1360 lsize = min(SIZE(x), SIZE(y))
1361 ALLOCATE(this%coord(lsize))
1362 this%coord(1:lsize)%x = x(1:lsize)
1363 this%coord(1:lsize)%y = y(1:lsize)
1364ENDIF
1365this%topo = optio_l(topo)
1367
1368END FUNCTION georef_coord_array_new
1369
1370
1371SUBROUTINE georef_coord_array_delete(this)
1372TYPE(georef_coord_array),INTENT(inout) :: this
1373
1374TYPE(georef_coord_array) :: lobj
1375
1376this = lobj
1377
1378END SUBROUTINE georef_coord_array_delete
1379
1380
1381ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1382TYPE(georef_coord_array),INTENT(in) :: this
1383LOGICAL :: res
1384
1385res = ALLOCATED(this%coord)
1386
1387END FUNCTION georef_coord_array_c_e
1388
1389
1394SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1395TYPE(georef_coord_array),INTENT(in) :: this
1396DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1397DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1398! allocatable per vedere di nascosto l'effetto che fa
1399INTEGER,OPTIONAL,INTENT(out) :: topo
1400TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1401
1402
1403IF (PRESENT(x)) THEN
1404 IF (ALLOCATED(this%coord)) THEN
1405 x = this%coord%x
1406 ENDIF
1407ENDIF
1408IF (PRESENT(y)) THEN
1409 IF (ALLOCATED(this%coord)) THEN
1410 y = this%coord%y
1411 ENDIF
1412ENDIF
1413IF (PRESENT(topo)) topo = this%topo
1415
1416END SUBROUTINE georef_coord_array_getval
1417
1418
1424SUBROUTINE georef_coord_array_compute_bbox(this)
1425TYPE(georef_coord_array),INTENT(inout) :: this
1426
1427IF (ALLOCATED(this%coord)) THEN
1428 this%bbox(1)%x = minval(this%coord(:)%x)
1429 this%bbox(1)%y = minval(this%coord(:)%y)
1430 this%bbox(2)%x = maxval(this%coord(:)%x)
1431 this%bbox(2)%y = maxval(this%coord(:)%y)
1432 this%bbox_updated = .true.
1433ENDIF
1434
1435END SUBROUTINE georef_coord_array_compute_bbox
1436
1437#ifdef HAVE_SHAPELIB
1438! internal method for importing a single shape
1439SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1440TYPE(georef_coord_array),INTENT(OUT) :: this
1441TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1442INTEGER,INTENT(IN) :: nshp
1443
1444TYPE(shpobject) :: shpobj
1445
1446! read shape object
1447shpobj = shpreadobject(shphandle, nshp)
1448IF (.NOT.shpisnull(shpobj)) THEN
1449! import it in georef_coord object
1450 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1451 topo=shpobj%nshptype)
1452 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1453 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1454 ELSE IF (ALLOCATED(this%parts)) THEN
1455 DEALLOCATE(this%parts)
1456 ENDIF
1457 CALL shpdestroyobject(shpobj)
1458 CALL compute_bbox(this)
1459ENDIF
1460
1461
1462END SUBROUTINE georef_coord_array_import
1463
1464
1465! internal method for exporting a single shape
1466SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1467TYPE(georef_coord_array),INTENT(in) :: this
1468TYPE(shpfileobject),INTENT(inout) :: shphandle
1469INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1470
1471INTEGER :: i
1472TYPE(shpobject) :: shpobj
1473
1474IF (ALLOCATED(this%coord)) THEN
1475 IF (ALLOCATED(this%parts)) THEN
1476 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1477 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1478 ELSE
1479 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1480 this%coord(:)%x, this%coord(:)%y)
1481 ENDIF
1482ELSE
1483 RETURN
1484ENDIF
1485
1486IF (.NOT.shpisnull(shpobj)) THEN
1487 i = shpwriteobject(shphandle, nshp, shpobj)
1488 CALL shpdestroyobject(shpobj)
1489ENDIF
1490
1491END SUBROUTINE georef_coord_array_export
1492
1493
1504SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1505TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1506CHARACTER(len=*),INTENT(in) :: shpfile
1507
1508REAL(kind=fp_d) :: minb(4), maxb(4)
1509INTEGER :: i, ns, shptype, dbfnf, dbfnr
1510TYPE(shpfileobject) :: shphandle
1511
1512shphandle = shpopen(trim(shpfile), 'rb')
1513IF (shpfileisnull(shphandle)) THEN
1514 ! log here
1515 CALL raise_error()
1516 RETURN
1517ENDIF
1518
1519! get info about file
1520CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1521IF (ns > 0) THEN ! allocate and read the object
1523 DO i = 1, ns
1524 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1525 ENDDO
1526ENDIF
1527
1528CALL shpclose(shphandle)
1529! pack object to save memory
1531
1532END SUBROUTINE arrayof_georef_coord_array_import
1533
1534
1540SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1541TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1542CHARACTER(len=*),INTENT(in) :: shpfile
1543
1544INTEGER :: i
1545TYPE(shpfileobject) :: shphandle
1546
1547IF (this%arraysize > 0) THEN
1548 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1549ELSE
1550 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1551ENDIF
1552IF (shpfileisnull(shphandle)) THEN
1553 ! log here
1554 CALL raise_error()
1555 RETURN
1556ENDIF
1557
1558DO i = 1, this%arraysize
1559 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1560ENDDO
1561
1562CALL shpclose(shphandle)
1563
1564END SUBROUTINE arrayof_georef_coord_array_export
1565#endif
1566
1578FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1579TYPE(georef_coord), INTENT(IN) :: this
1580TYPE(georef_coord_array), INTENT(IN) :: poly
1581LOGICAL :: inside
1582
1583INTEGER :: i
1584
1585inside = .false.
1587IF (.NOT.ALLOCATED(poly%coord)) RETURN
1588! if outside bounding box stop here
1589IF (poly%bbox_updated) THEN
1590 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1591ENDIF
1592
1593IF (ALLOCATED(poly%parts)) THEN
1594 DO i = 1, SIZE(poly%parts)-1
1596 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1597 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1598 ENDDO
1599 IF (SIZE(poly%parts) > 0) THEN ! safety check
1601 poly%coord(poly%parts(i)+1:)%x, &
1602 poly%coord(poly%parts(i)+1:)%y)
1603 ENDIF
1604
1605ELSE
1606 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1607 inside = pointinpoly(this%x, this%y, &
1608 poly%coord(:)%x, poly%coord(:)%y)
1609ENDIF
1610
1611CONTAINS
1612
1613FUNCTION pointinpoly(x, y, px, py)
1614DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1615LOGICAL :: pointinpoly
1616
1617INTEGER :: i, j, starti
1618
1619pointinpoly = .false.
1620
1621IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1622 starti = 2
1623 j = 1
1624ELSE ! unclosed polygon
1625 starti = 1
1626 j = SIZE(px)
1627ENDIF
1628DO i = starti, SIZE(px)
1629 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1630 (py(j) <= y .AND. y < py(i))) THEN
1631 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1632 pointinpoly = .NOT. pointinpoly
1633 ENDIF
1634 ENDIF
1635 j = i
1636ENDDO
1637
1638END FUNCTION pointinpoly
1639
1640END FUNCTION georef_coord_inside
1641
1642
1643
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 |