libsim Versione 7.1.11
|
◆ 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 927 del file georef_coord_class.F90. 928! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
929! authors:
930! Davide Cesari <dcesari@arpa.emr.it>
931! Paolo Patruno <ppatruno@arpa.emr.it>
932
933! This program is free software; you can redistribute it and/or
934! modify it under the terms of the GNU General Public License as
935! published by the Free Software Foundation; either version 2 of
936! the License, or (at your option) any later version.
937
938! This program is distributed in the hope that it will be useful,
939! but WITHOUT ANY WARRANTY; without even the implied warranty of
940! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
941! GNU General Public License for more details.
942
943! You should have received a copy of the GNU General Public License
944! along with this program. If not, see <http://www.gnu.org/licenses/>.
945#include "config.h"
946
964USE geo_proj_class
965#ifdef HAVE_SHAPELIB
966USE shapelib
967#endif
968IMPLICIT NONE
969
975 PRIVATE
976 DOUBLE PRECISION :: x=dmiss, y=dmiss
978
981
987 PRIVATE
988 INTEGER,ALLOCATABLE :: parts(:)
989 TYPE(georef_coord),ALLOCATABLE :: coord(:)
990 INTEGER :: topo=imiss
991 TYPE(geo_proj) :: proj
992 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
993 LOGICAL :: bbox_updated=.false.
995
996INTEGER,PARAMETER :: georef_coord_array_point = 1
997INTEGER,PARAMETER :: georef_coord_array_arc = 3
998INTEGER,PARAMETER :: georef_coord_array_polygon = 5
999INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
1000
1001
1006 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1007END INTERFACE
1008
1011 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1012END INTERFACE
1013
1016 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
1017END INTERFACE
1018
1019INTERFACE compute_bbox
1020 MODULE PROCEDURE georef_coord_array_compute_bbox
1021END INTERFACE
1022
1024INTERFACE OPERATOR (==)
1025 MODULE PROCEDURE georef_coord_eq
1026END INTERFACE
1027
1029INTERFACE OPERATOR (/=)
1030 MODULE PROCEDURE georef_coord_ne
1031END INTERFACE
1032
1035INTERFACE OPERATOR (>=)
1036 MODULE PROCEDURE georef_coord_ge
1037END INTERFACE
1038
1041INTERFACE OPERATOR (<=)
1042 MODULE PROCEDURE georef_coord_le
1043END INTERFACE
1044
1045#ifdef HAVE_SHAPELIB
1046
1048INTERFACE import
1049 MODULE PROCEDURE arrayof_georef_coord_array_import
1050END INTERFACE
1051
1055 MODULE PROCEDURE arrayof_georef_coord_array_export
1056END INTERFACE
1057#endif
1058
1062 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1063END INTERFACE
1064
1068 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1069END INTERFACE
1070
1073 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1074END INTERFACE
1075
1078 MODULE PROCEDURE georef_coord_dist
1079END INTERFACE
1080
1081#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1082#define ARRAYOF_TYPE arrayof_georef_coord_array
1083!define ARRAYOF_ORIGEQ 0
1084#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1085#include "arrayof_pre.F90"
1086! from arrayof
1088
1089PRIVATE
1091 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1092 georef_coord_array_polygon, georef_coord_array_multipoint, &
1094#ifdef HAVE_SHAPELIB
1096#endif
1098 georef_coord_new, georef_coord_array_new
1099
1100CONTAINS
1101
1102#include "arrayof_post.F90"
1103
1104! ===================
1105! == georef_coord ==
1106! ===================
1110FUNCTION georef_coord_new(x, y) RESULT(this)
1111DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1112DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1113TYPE(georef_coord) :: this
1114
1117
1118END FUNCTION georef_coord_new
1119
1120
1121SUBROUTINE georef_coord_delete(this)
1122TYPE(georef_coord),INTENT(inout) :: this
1123
1124this%x = dmiss
1125this%y = dmiss
1126
1127END SUBROUTINE georef_coord_delete
1128
1129
1130ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1131TYPE(georef_coord),INTENT(in) :: this
1132LOGICAL :: res
1133
1134res = .NOT. this == georef_coord_miss
1135
1136END FUNCTION georef_coord_c_e
1137
1138
1145ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1146TYPE(georef_coord),INTENT(in) :: this
1147DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1148DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1149
1150IF (PRESENT(x)) x = this%x
1151IF (PRESENT(y)) y = this%y
1152
1153END SUBROUTINE georef_coord_getval
1154
1155
1164ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1165TYPE(georef_coord),INTENT(in) :: this
1166TYPE(geo_proj),INTENT(in) :: proj
1167DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1168DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1169DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1170DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1171
1172DOUBLE PRECISION :: llon, llat
1173
1174IF (PRESENT(x)) x = this%x
1175IF (PRESENT(y)) y = this%y
1176IF (PRESENT(lon) .OR. present(lat)) THEN
1178 IF (PRESENT(lon)) lon = llon
1179 IF (PRESENT(lat)) lat = llat
1180ENDIF
1181
1182END SUBROUTINE georef_coord_proj_getval
1183
1184
1185! document and improve
1186ELEMENTAL FUNCTION getlat(this)
1187TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1188DOUBLE PRECISION :: getlat ! latitudine geografica
1189
1190getlat = this%y ! change!!!
1191
1192END FUNCTION getlat
1193
1194! document and improve
1195ELEMENTAL FUNCTION getlon(this)
1196TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1197DOUBLE PRECISION :: getlon ! longitudine geografica
1198
1199getlon = this%x ! change!!!
1200
1201END FUNCTION getlon
1202
1203
1204ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1205TYPE(georef_coord),INTENT(IN) :: this, that
1206LOGICAL :: res
1207
1208res = (this%x == that%x .AND. this%y == that%y)
1209
1210END FUNCTION georef_coord_eq
1211
1212
1213ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1214TYPE(georef_coord),INTENT(IN) :: this, that
1215LOGICAL :: res
1216
1217res = (this%x >= that%x .AND. this%y >= that%y)
1218
1219END FUNCTION georef_coord_ge
1220
1221
1222ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1223TYPE(georef_coord),INTENT(IN) :: this, that
1224LOGICAL :: res
1225
1226res = (this%x <= that%x .AND. this%y <= that%y)
1227
1228END FUNCTION georef_coord_le
1229
1230
1231ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1232TYPE(georef_coord),INTENT(IN) :: this, that
1233LOGICAL :: res
1234
1235res = .NOT.(this == that)
1236
1237END FUNCTION georef_coord_ne
1238
1239
1245SUBROUTINE georef_coord_read_unit(this, unit)
1246TYPE(georef_coord),INTENT(out) :: this
1247INTEGER, INTENT(in) :: unit
1248
1249CALL georef_coord_vect_read_unit((/this/), unit)
1250
1251END SUBROUTINE georef_coord_read_unit
1252
1253
1259SUBROUTINE georef_coord_vect_read_unit(this, unit)
1260TYPE(georef_coord) :: this(:)
1261INTEGER, INTENT(in) :: unit
1262
1263CHARACTER(len=40) :: form
1264INTEGER :: i
1265
1266INQUIRE(unit, form=form)
1267IF (form == 'FORMATTED') THEN
1268 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1269!TODO bug gfortran compiler !
1270!missing values are unredeable when formatted
1271ELSE
1272 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1273ENDIF
1274
1275END SUBROUTINE georef_coord_vect_read_unit
1276
1277
1282SUBROUTINE georef_coord_write_unit(this, unit)
1283TYPE(georef_coord),INTENT(in) :: this
1284INTEGER, INTENT(in) :: unit
1285
1286CALL georef_coord_vect_write_unit((/this/), unit)
1287
1288END SUBROUTINE georef_coord_write_unit
1289
1290
1295SUBROUTINE georef_coord_vect_write_unit(this, unit)
1296TYPE(georef_coord),INTENT(in) :: this(:)
1297INTEGER, INTENT(in) :: unit
1298
1299CHARACTER(len=40) :: form
1300INTEGER :: i
1301
1302INQUIRE(unit, form=form)
1303IF (form == 'FORMATTED') THEN
1304 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1305!TODO bug gfortran compiler !
1306!missing values are unredeable when formatted
1307ELSE
1308 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1309ENDIF
1310
1311END SUBROUTINE georef_coord_vect_write_unit
1312
1313
1316FUNCTION georef_coord_dist(this, that) RESULT(dist)
1318TYPE(georef_coord), INTENT (IN) :: this
1319TYPE(georef_coord), INTENT (IN) :: that
1320DOUBLE PRECISION :: dist
1321
1322DOUBLE PRECISION :: x,y
1323! Distanza approssimata, valida per piccoli angoli
1324
1325x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1326y = (this%y-that%y)
1327dist = sqrt(x**2 + y**2)*degrad*rearth
1328
1329END FUNCTION georef_coord_dist
1330
1331
1337FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1338TYPE(georef_coord),INTENT(IN) :: this
1339TYPE(georef_coord),INTENT(IN) :: coordmin
1340TYPE(georef_coord),INTENT(IN) :: coordmax
1341LOGICAL :: res
1342
1343res = (this >= coordmin .AND. this <= coordmax)
1344
1345END FUNCTION georef_coord_inside_rectang
1346
1347
1348! ========================
1349! == georef_coord_array ==
1350! ========================
1356FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1357DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1358DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1359INTEGER,INTENT(in),OPTIONAL :: topo
1360TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1361TYPE(georef_coord_array) :: this
1362
1363INTEGER :: lsize
1364
1365IF (PRESENT(x) .AND. PRESENT(y)) THEN
1366 lsize = min(SIZE(x), SIZE(y))
1367 ALLOCATE(this%coord(lsize))
1368 this%coord(1:lsize)%x = x(1:lsize)
1369 this%coord(1:lsize)%y = y(1:lsize)
1370ENDIF
1371this%topo = optio_l(topo)
1373
1374END FUNCTION georef_coord_array_new
1375
1376
1377SUBROUTINE georef_coord_array_delete(this)
1378TYPE(georef_coord_array),INTENT(inout) :: this
1379
1380TYPE(georef_coord_array) :: lobj
1381
1382this = lobj
1383
1384END SUBROUTINE georef_coord_array_delete
1385
1386
1387ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1388TYPE(georef_coord_array),INTENT(in) :: this
1389LOGICAL :: res
1390
1391res = ALLOCATED(this%coord)
1392
1393END FUNCTION georef_coord_array_c_e
1394
1395
1400SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1401TYPE(georef_coord_array),INTENT(in) :: this
1402DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1403DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1404! allocatable per vedere di nascosto l'effetto che fa
1405INTEGER,OPTIONAL,INTENT(out) :: topo
1406TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1407
1408
1409IF (PRESENT(x)) THEN
1410 IF (ALLOCATED(this%coord)) THEN
1411 x = this%coord%x
1412 ENDIF
1413ENDIF
1414IF (PRESENT(y)) THEN
1415 IF (ALLOCATED(this%coord)) THEN
1416 y = this%coord%y
1417 ENDIF
1418ENDIF
1419IF (PRESENT(topo)) topo = this%topo
1421
1422END SUBROUTINE georef_coord_array_getval
1423
1424
1430SUBROUTINE georef_coord_array_compute_bbox(this)
1431TYPE(georef_coord_array),INTENT(inout) :: this
1432
1433IF (ALLOCATED(this%coord)) THEN
1434 this%bbox(1)%x = minval(this%coord(:)%x)
1435 this%bbox(1)%y = minval(this%coord(:)%y)
1436 this%bbox(2)%x = maxval(this%coord(:)%x)
1437 this%bbox(2)%y = maxval(this%coord(:)%y)
1438 this%bbox_updated = .true.
1439ENDIF
1440
1441END SUBROUTINE georef_coord_array_compute_bbox
1442
1443#ifdef HAVE_SHAPELIB
1444! internal method for importing a single shape
1445SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1446TYPE(georef_coord_array),INTENT(OUT) :: this
1447TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1448INTEGER,INTENT(IN) :: nshp
1449
1450TYPE(shpobject) :: shpobj
1451
1452! read shape object
1453shpobj = shpreadobject(shphandle, nshp)
1454IF (.NOT.shpisnull(shpobj)) THEN
1455! import it in georef_coord object
1456 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1457 topo=shpobj%nshptype)
1458 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1459 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1460 ELSE IF (ALLOCATED(this%parts)) THEN
1461 DEALLOCATE(this%parts)
1462 ENDIF
1463 CALL shpdestroyobject(shpobj)
1464 CALL compute_bbox(this)
1465ENDIF
1466
1467
1468END SUBROUTINE georef_coord_array_import
1469
1470
1471! internal method for exporting a single shape
1472SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1473TYPE(georef_coord_array),INTENT(in) :: this
1474TYPE(shpfileobject),INTENT(inout) :: shphandle
1475INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1476
1477INTEGER :: i
1478TYPE(shpobject) :: shpobj
1479
1480IF (ALLOCATED(this%coord)) THEN
1481 IF (ALLOCATED(this%parts)) THEN
1482 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1483 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1484 ELSE
1485 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1486 this%coord(:)%x, this%coord(:)%y)
1487 ENDIF
1488ELSE
1489 RETURN
1490ENDIF
1491
1492IF (.NOT.shpisnull(shpobj)) THEN
1493 i = shpwriteobject(shphandle, nshp, shpobj)
1494 CALL shpdestroyobject(shpobj)
1495ENDIF
1496
1497END SUBROUTINE georef_coord_array_export
1498
1499
1510SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1511TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1512CHARACTER(len=*),INTENT(in) :: shpfile
1513
1514REAL(kind=fp_d) :: minb(4), maxb(4)
1515INTEGER :: i, ns, shptype, dbfnf, dbfnr
1516TYPE(shpfileobject) :: shphandle
1517
1518shphandle = shpopen(trim(shpfile), 'rb')
1519IF (shpfileisnull(shphandle)) THEN
1520 ! log here
1521 CALL raise_error()
1522 RETURN
1523ENDIF
1524
1525! get info about file
1526CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1527IF (ns > 0) THEN ! allocate and read the object
1529 DO i = 1, ns
1530 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1531 ENDDO
1532ENDIF
1533
1534CALL shpclose(shphandle)
1535! pack object to save memory
1537
1538END SUBROUTINE arrayof_georef_coord_array_import
1539
1540
1546SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1547TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1548CHARACTER(len=*),INTENT(in) :: shpfile
1549
1550INTEGER :: i
1551TYPE(shpfileobject) :: shphandle
1552
1553IF (this%arraysize > 0) THEN
1554 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1555ELSE
1556 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1557ENDIF
1558IF (shpfileisnull(shphandle)) THEN
1559 ! log here
1560 CALL raise_error()
1561 RETURN
1562ENDIF
1563
1564DO i = 1, this%arraysize
1565 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1566ENDDO
1567
1568CALL shpclose(shphandle)
1569
1570END SUBROUTINE arrayof_georef_coord_array_export
1571#endif
1572
1584FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1585TYPE(georef_coord), INTENT(IN) :: this
1586TYPE(georef_coord_array), INTENT(IN) :: poly
1587LOGICAL :: inside
1588
1589INTEGER :: i
1590
1591inside = .false.
1593IF (.NOT.ALLOCATED(poly%coord)) RETURN
1594! if outside bounding box stop here
1595IF (poly%bbox_updated) THEN
1596 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1597ENDIF
1598
1599IF (ALLOCATED(poly%parts)) THEN
1600 DO i = 1, SIZE(poly%parts)-1
1602 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1603 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1604 ENDDO
1605 IF (SIZE(poly%parts) > 0) THEN ! safety check
1607 poly%coord(poly%parts(i)+1:)%x, &
1608 poly%coord(poly%parts(i)+1:)%y)
1609 ENDIF
1610
1611ELSE
1612 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1613 inside = pointinpoly(this%x, this%y, &
1614 poly%coord(:)%x, poly%coord(:)%y)
1615ENDIF
1616
1617CONTAINS
1618
1619FUNCTION pointinpoly(x, y, px, py)
1620DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1621LOGICAL :: pointinpoly
1622
1623INTEGER :: i, j, starti
1624
1625pointinpoly = .false.
1626
1627IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1628 starti = 2
1629 j = 1
1630ELSE ! unclosed polygon
1631 starti = 1
1632 j = SIZE(px)
1633ENDIF
1634DO i = starti, SIZE(px)
1635 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1636 (py(j) <= y .AND. y < py(i))) THEN
1637 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1638 pointinpoly = .NOT. pointinpoly
1639 ENDIF
1640 ENDIF
1641 j = i
1642ENDDO
1643
1644END FUNCTION pointinpoly
1645
1646END FUNCTION georef_coord_inside
1647
1648
1649
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 |