libsim Versione 7.2.1
|
◆ 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 942 del file georef_coord_class.F90. 943! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
944! authors:
945! Davide Cesari <dcesari@arpa.emr.it>
946! Paolo Patruno <ppatruno@arpa.emr.it>
947
948! This program is free software; you can redistribute it and/or
949! modify it under the terms of the GNU General Public License as
950! published by the Free Software Foundation; either version 2 of
951! the License, or (at your option) any later version.
952
953! This program is distributed in the hope that it will be useful,
954! but WITHOUT ANY WARRANTY; without even the implied warranty of
955! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
956! GNU General Public License for more details.
957
958! You should have received a copy of the GNU General Public License
959! along with this program. If not, see <http://www.gnu.org/licenses/>.
960#include "config.h"
961
979USE geo_proj_class
980#ifdef HAVE_SHAPELIB
981USE shapelib
982#endif
983IMPLICIT NONE
984
990 PRIVATE
991 DOUBLE PRECISION :: x=dmiss, y=dmiss
993
996
1002 PRIVATE
1003 INTEGER,ALLOCATABLE :: parts(:)
1004 TYPE(georef_coord),ALLOCATABLE :: coord(:)
1005 INTEGER :: topo=imiss
1006 TYPE(geo_proj) :: proj
1007 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
1008 LOGICAL :: bbox_updated=.false.
1010
1011INTEGER,PARAMETER :: georef_coord_array_point = 1
1012INTEGER,PARAMETER :: georef_coord_array_arc = 3
1013INTEGER,PARAMETER :: georef_coord_array_polygon = 5
1014INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
1015
1016
1021 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1022END INTERFACE
1023
1026 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1027END INTERFACE
1028
1031 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
1032END INTERFACE
1033
1034INTERFACE compute_bbox
1035 MODULE PROCEDURE georef_coord_array_compute_bbox
1036END INTERFACE
1037
1039INTERFACE OPERATOR (==)
1040 MODULE PROCEDURE georef_coord_eq
1041END INTERFACE
1042
1044INTERFACE OPERATOR (/=)
1045 MODULE PROCEDURE georef_coord_ne
1046END INTERFACE
1047
1050INTERFACE OPERATOR (>=)
1051 MODULE PROCEDURE georef_coord_ge
1052END INTERFACE
1053
1056INTERFACE OPERATOR (<=)
1057 MODULE PROCEDURE georef_coord_le
1058END INTERFACE
1059
1060#ifdef HAVE_SHAPELIB
1061
1063INTERFACE import
1064 MODULE PROCEDURE arrayof_georef_coord_array_import
1065END INTERFACE
1066
1070 MODULE PROCEDURE arrayof_georef_coord_array_export
1071END INTERFACE
1072#endif
1073
1077 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1078END INTERFACE
1079
1083 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1084END INTERFACE
1085
1088 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1089END INTERFACE
1090
1093 MODULE PROCEDURE georef_coord_dist
1094END INTERFACE
1095
1096#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1097#define ARRAYOF_TYPE arrayof_georef_coord_array
1098!define ARRAYOF_ORIGEQ 0
1099#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1100#include "arrayof_pre.F90"
1101! from arrayof
1103
1104PRIVATE
1106 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1107 georef_coord_array_polygon, georef_coord_array_multipoint, &
1109#ifdef HAVE_SHAPELIB
1111#endif
1113 georef_coord_new, georef_coord_array_new
1114
1115CONTAINS
1116
1117#include "arrayof_post.F90"
1118
1119! ===================
1120! == georef_coord ==
1121! ===================
1125FUNCTION georef_coord_new(x, y) RESULT(this)
1126DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1127DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1128TYPE(georef_coord) :: this
1129
1132
1133END FUNCTION georef_coord_new
1134
1135
1136SUBROUTINE georef_coord_delete(this)
1137TYPE(georef_coord),INTENT(inout) :: this
1138
1139this%x = dmiss
1140this%y = dmiss
1141
1142END SUBROUTINE georef_coord_delete
1143
1144
1145ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1146TYPE(georef_coord),INTENT(in) :: this
1147LOGICAL :: res
1148
1149res = .NOT. this == georef_coord_miss
1150
1151END FUNCTION georef_coord_c_e
1152
1153
1160ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1161TYPE(georef_coord),INTENT(in) :: this
1162DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1163DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1164
1165IF (PRESENT(x)) x = this%x
1166IF (PRESENT(y)) y = this%y
1167
1168END SUBROUTINE georef_coord_getval
1169
1170
1179ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1180TYPE(georef_coord),INTENT(in) :: this
1181TYPE(geo_proj),INTENT(in) :: proj
1182DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1183DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1184DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1185DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1186
1187DOUBLE PRECISION :: llon, llat
1188
1189IF (PRESENT(x)) x = this%x
1190IF (PRESENT(y)) y = this%y
1191IF (PRESENT(lon) .OR. present(lat)) THEN
1193 IF (PRESENT(lon)) lon = llon
1194 IF (PRESENT(lat)) lat = llat
1195ENDIF
1196
1197END SUBROUTINE georef_coord_proj_getval
1198
1199
1200! document and improve
1201ELEMENTAL FUNCTION getlat(this)
1202TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1203DOUBLE PRECISION :: getlat ! latitudine geografica
1204
1205getlat = this%y ! change!!!
1206
1207END FUNCTION getlat
1208
1209! document and improve
1210ELEMENTAL FUNCTION getlon(this)
1211TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1212DOUBLE PRECISION :: getlon ! longitudine geografica
1213
1214getlon = this%x ! change!!!
1215
1216END FUNCTION getlon
1217
1218
1219ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1220TYPE(georef_coord),INTENT(IN) :: this, that
1221LOGICAL :: res
1222
1223res = (this%x == that%x .AND. this%y == that%y)
1224
1225END FUNCTION georef_coord_eq
1226
1227
1228ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1229TYPE(georef_coord),INTENT(IN) :: this, that
1230LOGICAL :: res
1231
1232res = (this%x >= that%x .AND. this%y >= that%y)
1233
1234END FUNCTION georef_coord_ge
1235
1236
1237ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1238TYPE(georef_coord),INTENT(IN) :: this, that
1239LOGICAL :: res
1240
1241res = (this%x <= that%x .AND. this%y <= that%y)
1242
1243END FUNCTION georef_coord_le
1244
1245
1246ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1247TYPE(georef_coord),INTENT(IN) :: this, that
1248LOGICAL :: res
1249
1250res = .NOT.(this == that)
1251
1252END FUNCTION georef_coord_ne
1253
1254
1260SUBROUTINE georef_coord_read_unit(this, unit)
1261TYPE(georef_coord),INTENT(out) :: this
1262INTEGER, INTENT(in) :: unit
1263
1264CALL georef_coord_vect_read_unit((/this/), unit)
1265
1266END SUBROUTINE georef_coord_read_unit
1267
1268
1274SUBROUTINE georef_coord_vect_read_unit(this, unit)
1275TYPE(georef_coord) :: this(:)
1276INTEGER, INTENT(in) :: unit
1277
1278CHARACTER(len=40) :: form
1279INTEGER :: i
1280
1281INQUIRE(unit, form=form)
1282IF (form == 'FORMATTED') THEN
1283 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1284!TODO bug gfortran compiler !
1285!missing values are unredeable when formatted
1286ELSE
1287 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1288ENDIF
1289
1290END SUBROUTINE georef_coord_vect_read_unit
1291
1292
1297SUBROUTINE georef_coord_write_unit(this, unit)
1298TYPE(georef_coord),INTENT(in) :: this
1299INTEGER, INTENT(in) :: unit
1300
1301CALL georef_coord_vect_write_unit((/this/), unit)
1302
1303END SUBROUTINE georef_coord_write_unit
1304
1305
1310SUBROUTINE georef_coord_vect_write_unit(this, unit)
1311TYPE(georef_coord),INTENT(in) :: this(:)
1312INTEGER, INTENT(in) :: unit
1313
1314CHARACTER(len=40) :: form
1315INTEGER :: i
1316
1317INQUIRE(unit, form=form)
1318IF (form == 'FORMATTED') THEN
1319 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1320!TODO bug gfortran compiler !
1321!missing values are unredeable when formatted
1322ELSE
1323 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1324ENDIF
1325
1326END SUBROUTINE georef_coord_vect_write_unit
1327
1328
1331FUNCTION georef_coord_dist(this, that) RESULT(dist)
1333TYPE(georef_coord), INTENT (IN) :: this
1334TYPE(georef_coord), INTENT (IN) :: that
1335DOUBLE PRECISION :: dist
1336
1337DOUBLE PRECISION :: x,y
1338! Distanza approssimata, valida per piccoli angoli
1339
1340x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1341y = (this%y-that%y)
1342dist = sqrt(x**2 + y**2)*degrad*rearth
1343
1344END FUNCTION georef_coord_dist
1345
1346
1352FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1353TYPE(georef_coord),INTENT(IN) :: this
1354TYPE(georef_coord),INTENT(IN) :: coordmin
1355TYPE(georef_coord),INTENT(IN) :: coordmax
1356LOGICAL :: res
1357
1358res = (this >= coordmin .AND. this <= coordmax)
1359
1360END FUNCTION georef_coord_inside_rectang
1361
1362
1363! ========================
1364! == georef_coord_array ==
1365! ========================
1371FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1372DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1373DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1374INTEGER,INTENT(in),OPTIONAL :: topo
1375TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1376TYPE(georef_coord_array) :: this
1377
1378INTEGER :: lsize
1379
1380IF (PRESENT(x) .AND. PRESENT(y)) THEN
1381 lsize = min(SIZE(x), SIZE(y))
1382 ALLOCATE(this%coord(lsize))
1383 this%coord(1:lsize)%x = x(1:lsize)
1384 this%coord(1:lsize)%y = y(1:lsize)
1385ENDIF
1386this%topo = optio_l(topo)
1388
1389END FUNCTION georef_coord_array_new
1390
1391
1392SUBROUTINE georef_coord_array_delete(this)
1393TYPE(georef_coord_array),INTENT(inout) :: this
1394
1395TYPE(georef_coord_array) :: lobj
1396
1397this = lobj
1398
1399END SUBROUTINE georef_coord_array_delete
1400
1401
1402ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1403TYPE(georef_coord_array),INTENT(in) :: this
1404LOGICAL :: res
1405
1406res = ALLOCATED(this%coord)
1407
1408END FUNCTION georef_coord_array_c_e
1409
1410
1415SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1416TYPE(georef_coord_array),INTENT(in) :: this
1417DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1418DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1419! allocatable per vedere di nascosto l'effetto che fa
1420INTEGER,OPTIONAL,INTENT(out) :: topo
1421TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1422
1423
1424IF (PRESENT(x)) THEN
1425 IF (ALLOCATED(this%coord)) THEN
1426 x = this%coord%x
1427 ENDIF
1428ENDIF
1429IF (PRESENT(y)) THEN
1430 IF (ALLOCATED(this%coord)) THEN
1431 y = this%coord%y
1432 ENDIF
1433ENDIF
1434IF (PRESENT(topo)) topo = this%topo
1436
1437END SUBROUTINE georef_coord_array_getval
1438
1439
1445SUBROUTINE georef_coord_array_compute_bbox(this)
1446TYPE(georef_coord_array),INTENT(inout) :: this
1447
1448IF (ALLOCATED(this%coord)) THEN
1449 this%bbox(1)%x = minval(this%coord(:)%x)
1450 this%bbox(1)%y = minval(this%coord(:)%y)
1451 this%bbox(2)%x = maxval(this%coord(:)%x)
1452 this%bbox(2)%y = maxval(this%coord(:)%y)
1453 this%bbox_updated = .true.
1454ENDIF
1455
1456END SUBROUTINE georef_coord_array_compute_bbox
1457
1458#ifdef HAVE_SHAPELIB
1459! internal method for importing a single shape
1460SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1461TYPE(georef_coord_array),INTENT(OUT) :: this
1462TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1463INTEGER,INTENT(IN) :: nshp
1464
1465TYPE(shpobject) :: shpobj
1466
1467! read shape object
1468shpobj = shpreadobject(shphandle, nshp)
1469IF (.NOT.shpisnull(shpobj)) THEN
1470! import it in georef_coord object
1471 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1472 topo=shpobj%nshptype)
1473 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1474 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1475 ELSE IF (ALLOCATED(this%parts)) THEN
1476 DEALLOCATE(this%parts)
1477 ENDIF
1478 CALL shpdestroyobject(shpobj)
1479 CALL compute_bbox(this)
1480ENDIF
1481
1482
1483END SUBROUTINE georef_coord_array_import
1484
1485
1486! internal method for exporting a single shape
1487SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1488TYPE(georef_coord_array),INTENT(in) :: this
1489TYPE(shpfileobject),INTENT(inout) :: shphandle
1490INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1491
1492INTEGER :: i
1493TYPE(shpobject) :: shpobj
1494
1495IF (ALLOCATED(this%coord)) THEN
1496 IF (ALLOCATED(this%parts)) THEN
1497 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1498 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1499 ELSE
1500 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1501 this%coord(:)%x, this%coord(:)%y)
1502 ENDIF
1503ELSE
1504 RETURN
1505ENDIF
1506
1507IF (.NOT.shpisnull(shpobj)) THEN
1508 i = shpwriteobject(shphandle, nshp, shpobj)
1509 CALL shpdestroyobject(shpobj)
1510ENDIF
1511
1512END SUBROUTINE georef_coord_array_export
1513
1514
1525SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1526TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1527CHARACTER(len=*),INTENT(in) :: shpfile
1528
1529REAL(kind=fp_d) :: minb(4), maxb(4)
1530INTEGER :: i, ns, shptype, dbfnf, dbfnr
1531TYPE(shpfileobject) :: shphandle
1532
1533shphandle = shpopen(trim(shpfile), 'rb')
1534IF (shpfileisnull(shphandle)) THEN
1535 ! log here
1536 CALL raise_error()
1537 RETURN
1538ENDIF
1539
1540! get info about file
1541CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1542IF (ns > 0) THEN ! allocate and read the object
1544 DO i = 1, ns
1545 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1546 ENDDO
1547ENDIF
1548
1549CALL shpclose(shphandle)
1550! pack object to save memory
1552
1553END SUBROUTINE arrayof_georef_coord_array_import
1554
1555
1561SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1562TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1563CHARACTER(len=*),INTENT(in) :: shpfile
1564
1565INTEGER :: i
1566TYPE(shpfileobject) :: shphandle
1567
1568IF (this%arraysize > 0) THEN
1569 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1570ELSE
1571 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1572ENDIF
1573IF (shpfileisnull(shphandle)) THEN
1574 ! log here
1575 CALL raise_error()
1576 RETURN
1577ENDIF
1578
1579DO i = 1, this%arraysize
1580 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1581ENDDO
1582
1583CALL shpclose(shphandle)
1584
1585END SUBROUTINE arrayof_georef_coord_array_export
1586#endif
1587
1599FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1600TYPE(georef_coord), INTENT(IN) :: this
1601TYPE(georef_coord_array), INTENT(IN) :: poly
1602LOGICAL :: inside
1603
1604INTEGER :: i
1605
1606inside = .false.
1608IF (.NOT.ALLOCATED(poly%coord)) RETURN
1609! if outside bounding box stop here
1610IF (poly%bbox_updated) THEN
1611 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1612ENDIF
1613
1614IF (ALLOCATED(poly%parts)) THEN
1615 DO i = 1, SIZE(poly%parts)-1
1617 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1618 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1619 ENDDO
1620 IF (SIZE(poly%parts) > 0) THEN ! safety check
1622 poly%coord(poly%parts(i)+1:)%x, &
1623 poly%coord(poly%parts(i)+1:)%y)
1624 ENDIF
1625
1626ELSE
1627 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1628 inside = pointinpoly(this%x, this%y, &
1629 poly%coord(:)%x, poly%coord(:)%y)
1630ENDIF
1631
1632CONTAINS
1633
1634FUNCTION pointinpoly(x, y, px, py)
1635DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1636LOGICAL :: pointinpoly
1637
1638INTEGER :: i, j, starti
1639
1640pointinpoly = .false.
1641
1642IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1643 starti = 2
1644 j = 1
1645ELSE ! unclosed polygon
1646 starti = 1
1647 j = SIZE(px)
1648ENDIF
1649DO i = starti, SIZE(px)
1650 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1651 (py(j) <= y .AND. y < py(i))) THEN
1652 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1653 pointinpoly = .NOT. pointinpoly
1654 ENDIF
1655 ENDIF
1656 j = i
1657ENDDO
1658
1659END FUNCTION pointinpoly
1660
1661END FUNCTION georef_coord_inside
1662
1663
1664
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 |