libsim Versione 7.1.11
|
◆ georef_coord_inside_rectang()
Determines whether the point this lies inside a specified rectangle. The rectangle is oriented parallely to the coordinate system, its lower-left and upper-right vertices are specified by the two arguments. The function returns
Definizione alla linea 948 del file georef_coord_class.F90. 949! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
950! authors:
951! Davide Cesari <dcesari@arpa.emr.it>
952! Paolo Patruno <ppatruno@arpa.emr.it>
953
954! This program is free software; you can redistribute it and/or
955! modify it under the terms of the GNU General Public License as
956! published by the Free Software Foundation; either version 2 of
957! the License, or (at your option) any later version.
958
959! This program is distributed in the hope that it will be useful,
960! but WITHOUT ANY WARRANTY; without even the implied warranty of
961! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
962! GNU General Public License for more details.
963
964! You should have received a copy of the GNU General Public License
965! along with this program. If not, see <http://www.gnu.org/licenses/>.
966#include "config.h"
967
985USE geo_proj_class
986#ifdef HAVE_SHAPELIB
987USE shapelib
988#endif
989IMPLICIT NONE
990
996 PRIVATE
997 DOUBLE PRECISION :: x=dmiss, y=dmiss
999
1002
1008 PRIVATE
1009 INTEGER,ALLOCATABLE :: parts(:)
1010 TYPE(georef_coord),ALLOCATABLE :: coord(:)
1011 INTEGER :: topo=imiss
1012 TYPE(geo_proj) :: proj
1013 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
1014 LOGICAL :: bbox_updated=.false.
1016
1017INTEGER,PARAMETER :: georef_coord_array_point = 1
1018INTEGER,PARAMETER :: georef_coord_array_arc = 3
1019INTEGER,PARAMETER :: georef_coord_array_polygon = 5
1020INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
1021
1022
1027 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1028END INTERFACE
1029
1032 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1033END INTERFACE
1034
1037 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
1038END INTERFACE
1039
1040INTERFACE compute_bbox
1041 MODULE PROCEDURE georef_coord_array_compute_bbox
1042END INTERFACE
1043
1045INTERFACE OPERATOR (==)
1046 MODULE PROCEDURE georef_coord_eq
1047END INTERFACE
1048
1050INTERFACE OPERATOR (/=)
1051 MODULE PROCEDURE georef_coord_ne
1052END INTERFACE
1053
1056INTERFACE OPERATOR (>=)
1057 MODULE PROCEDURE georef_coord_ge
1058END INTERFACE
1059
1062INTERFACE OPERATOR (<=)
1063 MODULE PROCEDURE georef_coord_le
1064END INTERFACE
1065
1066#ifdef HAVE_SHAPELIB
1067
1069INTERFACE import
1070 MODULE PROCEDURE arrayof_georef_coord_array_import
1071END INTERFACE
1072
1076 MODULE PROCEDURE arrayof_georef_coord_array_export
1077END INTERFACE
1078#endif
1079
1083 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1084END INTERFACE
1085
1089 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1090END INTERFACE
1091
1094 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1095END INTERFACE
1096
1099 MODULE PROCEDURE georef_coord_dist
1100END INTERFACE
1101
1102#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1103#define ARRAYOF_TYPE arrayof_georef_coord_array
1104!define ARRAYOF_ORIGEQ 0
1105#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1106#include "arrayof_pre.F90"
1107! from arrayof
1109
1110PRIVATE
1112 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1113 georef_coord_array_polygon, georef_coord_array_multipoint, &
1115#ifdef HAVE_SHAPELIB
1117#endif
1119 georef_coord_new, georef_coord_array_new
1120
1121CONTAINS
1122
1123#include "arrayof_post.F90"
1124
1125! ===================
1126! == georef_coord ==
1127! ===================
1131FUNCTION georef_coord_new(x, y) RESULT(this)
1132DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1133DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1134TYPE(georef_coord) :: this
1135
1138
1139END FUNCTION georef_coord_new
1140
1141
1142SUBROUTINE georef_coord_delete(this)
1143TYPE(georef_coord),INTENT(inout) :: this
1144
1145this%x = dmiss
1146this%y = dmiss
1147
1148END SUBROUTINE georef_coord_delete
1149
1150
1151ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1152TYPE(georef_coord),INTENT(in) :: this
1153LOGICAL :: res
1154
1155res = .NOT. this == georef_coord_miss
1156
1157END FUNCTION georef_coord_c_e
1158
1159
1166ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1167TYPE(georef_coord),INTENT(in) :: this
1168DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1169DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1170
1171IF (PRESENT(x)) x = this%x
1172IF (PRESENT(y)) y = this%y
1173
1174END SUBROUTINE georef_coord_getval
1175
1176
1185ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1186TYPE(georef_coord),INTENT(in) :: this
1187TYPE(geo_proj),INTENT(in) :: proj
1188DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1189DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1190DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1191DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1192
1193DOUBLE PRECISION :: llon, llat
1194
1195IF (PRESENT(x)) x = this%x
1196IF (PRESENT(y)) y = this%y
1197IF (PRESENT(lon) .OR. present(lat)) THEN
1199 IF (PRESENT(lon)) lon = llon
1200 IF (PRESENT(lat)) lat = llat
1201ENDIF
1202
1203END SUBROUTINE georef_coord_proj_getval
1204
1205
1206! document and improve
1207ELEMENTAL FUNCTION getlat(this)
1208TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1209DOUBLE PRECISION :: getlat ! latitudine geografica
1210
1211getlat = this%y ! change!!!
1212
1213END FUNCTION getlat
1214
1215! document and improve
1216ELEMENTAL FUNCTION getlon(this)
1217TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1218DOUBLE PRECISION :: getlon ! longitudine geografica
1219
1220getlon = this%x ! change!!!
1221
1222END FUNCTION getlon
1223
1224
1225ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1226TYPE(georef_coord),INTENT(IN) :: this, that
1227LOGICAL :: res
1228
1229res = (this%x == that%x .AND. this%y == that%y)
1230
1231END FUNCTION georef_coord_eq
1232
1233
1234ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1235TYPE(georef_coord),INTENT(IN) :: this, that
1236LOGICAL :: res
1237
1238res = (this%x >= that%x .AND. this%y >= that%y)
1239
1240END FUNCTION georef_coord_ge
1241
1242
1243ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1244TYPE(georef_coord),INTENT(IN) :: this, that
1245LOGICAL :: res
1246
1247res = (this%x <= that%x .AND. this%y <= that%y)
1248
1249END FUNCTION georef_coord_le
1250
1251
1252ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1253TYPE(georef_coord),INTENT(IN) :: this, that
1254LOGICAL :: res
1255
1256res = .NOT.(this == that)
1257
1258END FUNCTION georef_coord_ne
1259
1260
1266SUBROUTINE georef_coord_read_unit(this, unit)
1267TYPE(georef_coord),INTENT(out) :: this
1268INTEGER, INTENT(in) :: unit
1269
1270CALL georef_coord_vect_read_unit((/this/), unit)
1271
1272END SUBROUTINE georef_coord_read_unit
1273
1274
1280SUBROUTINE georef_coord_vect_read_unit(this, unit)
1281TYPE(georef_coord) :: this(:)
1282INTEGER, INTENT(in) :: unit
1283
1284CHARACTER(len=40) :: form
1285INTEGER :: i
1286
1287INQUIRE(unit, form=form)
1288IF (form == 'FORMATTED') THEN
1289 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1290!TODO bug gfortran compiler !
1291!missing values are unredeable when formatted
1292ELSE
1293 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1294ENDIF
1295
1296END SUBROUTINE georef_coord_vect_read_unit
1297
1298
1303SUBROUTINE georef_coord_write_unit(this, unit)
1304TYPE(georef_coord),INTENT(in) :: this
1305INTEGER, INTENT(in) :: unit
1306
1307CALL georef_coord_vect_write_unit((/this/), unit)
1308
1309END SUBROUTINE georef_coord_write_unit
1310
1311
1316SUBROUTINE georef_coord_vect_write_unit(this, unit)
1317TYPE(georef_coord),INTENT(in) :: this(:)
1318INTEGER, INTENT(in) :: unit
1319
1320CHARACTER(len=40) :: form
1321INTEGER :: i
1322
1323INQUIRE(unit, form=form)
1324IF (form == 'FORMATTED') THEN
1325 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1326!TODO bug gfortran compiler !
1327!missing values are unredeable when formatted
1328ELSE
1329 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1330ENDIF
1331
1332END SUBROUTINE georef_coord_vect_write_unit
1333
1334
1337FUNCTION georef_coord_dist(this, that) RESULT(dist)
1339TYPE(georef_coord), INTENT (IN) :: this
1340TYPE(georef_coord), INTENT (IN) :: that
1341DOUBLE PRECISION :: dist
1342
1343DOUBLE PRECISION :: x,y
1344! Distanza approssimata, valida per piccoli angoli
1345
1346x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1347y = (this%y-that%y)
1348dist = sqrt(x**2 + y**2)*degrad*rearth
1349
1350END FUNCTION georef_coord_dist
1351
1352
1358FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1359TYPE(georef_coord),INTENT(IN) :: this
1360TYPE(georef_coord),INTENT(IN) :: coordmin
1361TYPE(georef_coord),INTENT(IN) :: coordmax
1362LOGICAL :: res
1363
1364res = (this >= coordmin .AND. this <= coordmax)
1365
1366END FUNCTION georef_coord_inside_rectang
1367
1368
1369! ========================
1370! == georef_coord_array ==
1371! ========================
1377FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1378DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1379DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1380INTEGER,INTENT(in),OPTIONAL :: topo
1381TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1382TYPE(georef_coord_array) :: this
1383
1384INTEGER :: lsize
1385
1386IF (PRESENT(x) .AND. PRESENT(y)) THEN
1387 lsize = min(SIZE(x), SIZE(y))
1388 ALLOCATE(this%coord(lsize))
1389 this%coord(1:lsize)%x = x(1:lsize)
1390 this%coord(1:lsize)%y = y(1:lsize)
1391ENDIF
1392this%topo = optio_l(topo)
1394
1395END FUNCTION georef_coord_array_new
1396
1397
1398SUBROUTINE georef_coord_array_delete(this)
1399TYPE(georef_coord_array),INTENT(inout) :: this
1400
1401TYPE(georef_coord_array) :: lobj
1402
1403this = lobj
1404
1405END SUBROUTINE georef_coord_array_delete
1406
1407
1408ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1409TYPE(georef_coord_array),INTENT(in) :: this
1410LOGICAL :: res
1411
1412res = ALLOCATED(this%coord)
1413
1414END FUNCTION georef_coord_array_c_e
1415
1416
1421SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1422TYPE(georef_coord_array),INTENT(in) :: this
1423DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1424DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1425! allocatable per vedere di nascosto l'effetto che fa
1426INTEGER,OPTIONAL,INTENT(out) :: topo
1427TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1428
1429
1430IF (PRESENT(x)) THEN
1431 IF (ALLOCATED(this%coord)) THEN
1432 x = this%coord%x
1433 ENDIF
1434ENDIF
1435IF (PRESENT(y)) THEN
1436 IF (ALLOCATED(this%coord)) THEN
1437 y = this%coord%y
1438 ENDIF
1439ENDIF
1440IF (PRESENT(topo)) topo = this%topo
1442
1443END SUBROUTINE georef_coord_array_getval
1444
1445
1451SUBROUTINE georef_coord_array_compute_bbox(this)
1452TYPE(georef_coord_array),INTENT(inout) :: this
1453
1454IF (ALLOCATED(this%coord)) THEN
1455 this%bbox(1)%x = minval(this%coord(:)%x)
1456 this%bbox(1)%y = minval(this%coord(:)%y)
1457 this%bbox(2)%x = maxval(this%coord(:)%x)
1458 this%bbox(2)%y = maxval(this%coord(:)%y)
1459 this%bbox_updated = .true.
1460ENDIF
1461
1462END SUBROUTINE georef_coord_array_compute_bbox
1463
1464#ifdef HAVE_SHAPELIB
1465! internal method for importing a single shape
1466SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1467TYPE(georef_coord_array),INTENT(OUT) :: this
1468TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1469INTEGER,INTENT(IN) :: nshp
1470
1471TYPE(shpobject) :: shpobj
1472
1473! read shape object
1474shpobj = shpreadobject(shphandle, nshp)
1475IF (.NOT.shpisnull(shpobj)) THEN
1476! import it in georef_coord object
1477 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1478 topo=shpobj%nshptype)
1479 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1480 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1481 ELSE IF (ALLOCATED(this%parts)) THEN
1482 DEALLOCATE(this%parts)
1483 ENDIF
1484 CALL shpdestroyobject(shpobj)
1485 CALL compute_bbox(this)
1486ENDIF
1487
1488
1489END SUBROUTINE georef_coord_array_import
1490
1491
1492! internal method for exporting a single shape
1493SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1494TYPE(georef_coord_array),INTENT(in) :: this
1495TYPE(shpfileobject),INTENT(inout) :: shphandle
1496INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1497
1498INTEGER :: i
1499TYPE(shpobject) :: shpobj
1500
1501IF (ALLOCATED(this%coord)) THEN
1502 IF (ALLOCATED(this%parts)) THEN
1503 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1504 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1505 ELSE
1506 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1507 this%coord(:)%x, this%coord(:)%y)
1508 ENDIF
1509ELSE
1510 RETURN
1511ENDIF
1512
1513IF (.NOT.shpisnull(shpobj)) THEN
1514 i = shpwriteobject(shphandle, nshp, shpobj)
1515 CALL shpdestroyobject(shpobj)
1516ENDIF
1517
1518END SUBROUTINE georef_coord_array_export
1519
1520
1531SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1532TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1533CHARACTER(len=*),INTENT(in) :: shpfile
1534
1535REAL(kind=fp_d) :: minb(4), maxb(4)
1536INTEGER :: i, ns, shptype, dbfnf, dbfnr
1537TYPE(shpfileobject) :: shphandle
1538
1539shphandle = shpopen(trim(shpfile), 'rb')
1540IF (shpfileisnull(shphandle)) THEN
1541 ! log here
1542 CALL raise_error()
1543 RETURN
1544ENDIF
1545
1546! get info about file
1547CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1548IF (ns > 0) THEN ! allocate and read the object
1550 DO i = 1, ns
1551 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1552 ENDDO
1553ENDIF
1554
1555CALL shpclose(shphandle)
1556! pack object to save memory
1558
1559END SUBROUTINE arrayof_georef_coord_array_import
1560
1561
1567SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1568TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1569CHARACTER(len=*),INTENT(in) :: shpfile
1570
1571INTEGER :: i
1572TYPE(shpfileobject) :: shphandle
1573
1574IF (this%arraysize > 0) THEN
1575 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1576ELSE
1577 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1578ENDIF
1579IF (shpfileisnull(shphandle)) THEN
1580 ! log here
1581 CALL raise_error()
1582 RETURN
1583ENDIF
1584
1585DO i = 1, this%arraysize
1586 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1587ENDDO
1588
1589CALL shpclose(shphandle)
1590
1591END SUBROUTINE arrayof_georef_coord_array_export
1592#endif
1593
1605FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1606TYPE(georef_coord), INTENT(IN) :: this
1607TYPE(georef_coord_array), INTENT(IN) :: poly
1608LOGICAL :: inside
1609
1610INTEGER :: i
1611
1612inside = .false.
1614IF (.NOT.ALLOCATED(poly%coord)) RETURN
1615! if outside bounding box stop here
1616IF (poly%bbox_updated) THEN
1617 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1618ENDIF
1619
1620IF (ALLOCATED(poly%parts)) THEN
1621 DO i = 1, SIZE(poly%parts)-1
1623 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1624 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1625 ENDDO
1626 IF (SIZE(poly%parts) > 0) THEN ! safety check
1628 poly%coord(poly%parts(i)+1:)%x, &
1629 poly%coord(poly%parts(i)+1:)%y)
1630 ENDIF
1631
1632ELSE
1633 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1634 inside = pointinpoly(this%x, this%y, &
1635 poly%coord(:)%x, poly%coord(:)%y)
1636ENDIF
1637
1638CONTAINS
1639
1640FUNCTION pointinpoly(x, y, px, py)
1641DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1642LOGICAL :: pointinpoly
1643
1644INTEGER :: i, j, starti
1645
1646pointinpoly = .false.
1647
1648IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1649 starti = 2
1650 j = 1
1651ELSE ! unclosed polygon
1652 starti = 1
1653 j = SIZE(px)
1654ENDIF
1655DO i = starti, SIZE(px)
1656 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1657 (py(j) <= y .AND. y < py(i))) THEN
1658 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1659 pointinpoly = .NOT. pointinpoly
1660 ENDIF
1661 ENDIF
1662 j = i
1663ENDDO
1664
1665END FUNCTION pointinpoly
1666
1667END FUNCTION georef_coord_inside
1668
1669
1670
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 |