libsim Versione 7.2.6
|
◆ 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!> This module defines objects describing georeferenced sparse points
981!! possibly with topology and projection information. This module
982!! defines two classes, \a georef_coord, which represents a single
983!! georeferenced point on the Earth, and \a georef_coord_array which
984!! defines a set of points with a topological relation.
985!!
986!! Both classes have \a PRIVATE members, so that they cannot be
987!! manipulated directly, but only through the proper methods.
988!!
989!! It is also possible to dafine a dynamically extendible array of \a
990!! georef_coord_array objects, of type \a arrayof_georef_coord_array,
991!! suitable for importing/exporting data from/to a shapefile.
992!!
993!! \ingroup base
998USE geo_proj_class
999#ifdef HAVE_SHAPELIB
1000USE shapelib
1001#endif
1002IMPLICIT NONE
1003
1004!> Derive type defining a single georeferenced point, either in
1005!! geodetic or in projected coordinates. The object has no information
1006!! about the georeferenced coordinate system associated, it must be
1007!! kept separately by the user.
1009 PRIVATE
1010 DOUBLE PRECISION :: x=dmiss, y=dmiss
1012
1013!> Missing value for georef_coord
1015
1016!> Derived type defining a one-dimensional array of georeferenced points
1017!! with an associated topology (isolated point, arc, polygon, group of
1018!! points), possibly broken into parts and with an associated
1019!! georeferenced coordinate system.
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 !< Topology for georef_coord_array (from shapelib): isolated point
1031INTEGER,PARAMETER :: georef_coord_array_arc = 3 !< Topology for georef_coord_array (from shapelib): arc (multiple arcs unsupported)
1032INTEGER,PARAMETER :: georef_coord_array_polygon = 5 !< Topology for georef_coord_array (from shapelib): polygon (necessarily closed, multiple polygons unsupported)
1033INTEGER,PARAMETER :: georef_coord_array_multipoint = 8 !< Topology for georef_coord_array (from shapelib): group of points
1034
1035
1036!> Detructors for the two classes.
1037!! They clean up all the information associated with the corresponding
1038!! objects.
1040 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1041END INTERFACE
1042
1043!> Check missing value.
1045 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1046END INTERFACE
1047
1048!> Methods for returning the value of object members.
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
1057!> Logical equality operator.
1058INTERFACE OPERATOR (==)
1059 MODULE PROCEDURE georef_coord_eq
1060END INTERFACE
1061
1062!> Logical inequality operator.
1063INTERFACE OPERATOR (/=)
1064 MODULE PROCEDURE georef_coord_ne
1065END INTERFACE
1066
1067!> Logical greater-equal operator. Returns true if the first point
1068!! lies to the northeast of the second
1069INTERFACE OPERATOR (>=)
1070 MODULE PROCEDURE georef_coord_ge
1071END INTERFACE
1072
1073!> Logical less-equal operator. Returns true if the first point
1074!! lies to the southwest of the second
1075INTERFACE OPERATOR (<=)
1076 MODULE PROCEDURE georef_coord_le
1077END INTERFACE
1078
1079#ifdef HAVE_SHAPELIB
1080!> Import an array of \a georef_coord_array objects from a file
1081!! in ESRI/Shapefile format.
1082INTERFACE import
1083 MODULE PROCEDURE arrayof_georef_coord_array_import
1084END INTERFACE
1085
1086!> Export an array of \a georef_coord_array objects to a file
1087!! in ESRI/Shapefile format.
1089 MODULE PROCEDURE arrayof_georef_coord_array_export
1090END INTERFACE
1091#endif
1092
1093!> Read a single \a georef_coord object or an array of \a georef_coord objects
1094!! from a Fortran \c FORMATTED or \c UNFORMATTED file.
1096 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1097END INTERFACE
1098
1099!> Write a single \a georef_coord object or an array of \a georef_coord objects
1100!! to a Fortran \c FORMATTED or \c UNFORMATTED file.
1102 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1103END INTERFACE
1104
1105!> Determine whether a point lies inside a polygon or a rectangle.
1107 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1108END INTERFACE
1109
1110!> Compute the distance in m between two points
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! ===================
1141!> Construct a \a georef_coord object with the optional parameters provided.
1142!! If coordinates are not provided the object obtained is empty
1143!! (missing, see c_e function).
1144FUNCTION georef_coord_new(x, y) RESULT(this)
1145DOUBLE PRECISION,INTENT(in),OPTIONAL :: x !< x coordinate
1146DOUBLE PRECISION,INTENT(in),OPTIONAL :: y !< y coordinate
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
1173!> Query a \a georef_coord object.
1174!! This is the correct way to retrieve the contents of a \a
1175!! georef_coord object, since its members are declared as \a
1176!! PRIVATE. It is declared as \a ELEMENTAL, thus it works also on
1177!! arrays of any shape and, in that case, the result will hae the same
1178!! shape as \a this.
1179ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1180TYPE(georef_coord),INTENT(in) :: this !< object to query
1181DOUBLE PRECISION,INTENT(out),OPTIONAL :: x !< x-coordinate
1182DOUBLE PRECISION,INTENT(out),OPTIONAL :: y !< y-coordinate
1183
1184IF (PRESENT(x)) x = this%x
1185IF (PRESENT(y)) y = this%y
1186
1187END SUBROUTINE georef_coord_getval
1188
1189
1190!> Query a \a georef_coord object associating a geographical projection to it.
1191!! This method allow to interpret the \a x,y coordinates of a \a
1192!! georef_coord object as projected on a specified geographical
1193!! projection system and retrieve the geodetic longitude and/or
1194!! latitude associated to them. When \a x or \y are requested it works
1195!! as the basic get_val method. It is declared as \a ELEMENTAL, thus
1196!! it works also on arrays of any shape and, in that case, the result
1197!! will hae the same shape as \a this.
1198ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1199TYPE(georef_coord),INTENT(in) :: this !< object to query
1200TYPE(geo_proj),INTENT(in) :: proj !< geographical projection associated to coordinates of \a this
1201DOUBLE PRECISION,INTENT(out),OPTIONAL :: x !< x-coordinate
1202DOUBLE PRECISION,INTENT(out),OPTIONAL :: y !< y-coordinate
1203DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon !< geodetic longitude
1204DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat !< geodetic latitude
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
1274!> Legge da un'unità di file il contenuto dell'oggetto \a this.
1275!! Il record da leggere deve essere stato scritto con la ::write_unit
1276!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
1277!! del vettore devono essere accordate. Il metodo controlla se il file è
1278!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1279SUBROUTINE georef_coord_read_unit(this, unit)
1280TYPE(georef_coord),INTENT(out) :: this !< oggetto da leggere
1281INTEGER, INTENT(in) :: unit !< unità da cui leggere
1282
1283CALL georef_coord_vect_read_unit((/this/), unit)
1284
1285END SUBROUTINE georef_coord_read_unit
1286
1287
1288!> Legge da un'unità di file il contenuto dell'oggetto \a this.
1289!! Il record da leggere deve essere stato scritto con la ::write_unit
1290!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
1291!! del vettore devono essere accordate. Il metodo controlla se il file è
1292!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1293SUBROUTINE georef_coord_vect_read_unit(this, unit)
1294TYPE(georef_coord) :: this(:) !< oggetto da leggere
1295INTEGER, INTENT(in) :: unit !< unità da cui leggere
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
1312!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
1313!! Il record scritto potrà successivamente essere letto con la ::read_unit.
1314!! Il metodo controlla se il file è
1315!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1316SUBROUTINE georef_coord_write_unit(this, unit)
1317TYPE(georef_coord),INTENT(in) :: this !< oggetto da scrivere
1318INTEGER, INTENT(in) :: unit !< unità su cui scrivere
1319
1320CALL georef_coord_vect_write_unit((/this/), unit)
1321
1322END SUBROUTINE georef_coord_write_unit
1323
1324
1325!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
1326!! Il record scritto potrà successivamente essere letto con la ::read_unit.
1327!! Il metodo controlla se il file è
1328!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1329SUBROUTINE georef_coord_vect_write_unit(this, unit)
1330TYPE(georef_coord),INTENT(in) :: this(:) !< oggetto da scrivere
1331INTEGER, INTENT(in) :: unit !< unità su cui scrivere
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
1348!> Restituisce la distanza in m tra 2 oggetti georef_coord.
1349!! La distanza è calcolata approssimativamente ed è valida per piccoli angoli.
1350FUNCTION georef_coord_dist(this, that) RESULT(dist)
1352TYPE(georef_coord), INTENT (IN) :: this !< primo punto
1353TYPE(georef_coord), INTENT (IN) :: that !< secondo punto
1354DOUBLE PRECISION :: dist !< distanza in metri
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
1366!> Determines whether the point \a this lies inside a specified rectangle.
1367!! The rectangle is oriented parallely to the coordinate system, its
1368!! lower-left and upper-right vertices are specified by the two
1369!! arguments. The function returns \c .TRUE. also in the case it lies
1370!! exactly on the border of the rectangle.
1371FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1372TYPE(georef_coord),INTENT(IN) :: this !< point to test
1373TYPE(georef_coord),INTENT(IN) :: coordmin !< lower-left vertex
1374TYPE(georef_coord),INTENT(IN) :: coordmax !< upper-right vertex
1375LOGICAL :: res
1376
1377res = (this >= coordmin .AND. this <= coordmax)
1378
1379END FUNCTION georef_coord_inside_rectang
1380
1381
1382! ========================
1383! == georef_coord_array ==
1384! ========================
1385!> Construct a \a georef_coord_array object with the optional parameters provided.
1386!! If coordinates are not provided the object obtained is empty
1387!! (missing, see c_e function), if coordinate arrays are of different
1388!! lengths the \a georef_coord_array is initialised to the shortest
1389!! length provided.
1390FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1391DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:) !< x coordinate array
1392DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:) !< y coordinate array
1393INTEGER,INTENT(in),OPTIONAL :: topo !< topology of the object, one of the \a georef_coord_array_* constants
1394TYPE(geo_proj),INTENT(in),OPTIONAL :: proj !< geographical projection associated to the coordinates
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
1430!> Query a \a georef_coord_array object.
1431!! This is the correct way to retrieve the contents of a \a
1432!! georef_coord_array object, since its members are declared as \a
1433!! PRIVATE.
1434SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1435TYPE(georef_coord_array),INTENT(in) :: this !< object to query
1436DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:) !< x-coordinate
1437DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:) !< y-coordinate
1438! allocatable per vedere di nascosto l'effetto che fa
1439INTEGER,OPTIONAL,INTENT(out) :: topo !< topology associated with the coordinates
1440TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj !< geographical projection
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
1459!> Compute the bounding box of each shape in \a georef_coord_array object.
1460!! The bounding box is computed and stored in the object, it is used
1461!! by the inside() function for speedup; after it is computed the
1462!! object cannot be changed, otherwise the bounding box will not be
1463!! valid.
1464SUBROUTINE georef_coord_array_compute_bbox(this)
1465TYPE(georef_coord_array),INTENT(inout) :: this !< object to manipulate
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
1534!> Import an array of \a georef_coord_array objects from a file
1535!! in ESRI/Shapefile format. The \a this argument is an uninitialised
1536!! \a arrayof_georef_coord_array, every element of which,
1537!! this%array(n), is of type \a georef_coord_array and, on return,
1538!! will contain information from the n-th shape of the file. Topology
1539!! information and possible polygon parts are imported as well, while
1540!! no projection information, even if available, is imported. An
1541!! error condition while opening the file can be detected by checking
1542!! .NOT.ASSOCIATED(this%array), while an error reading shape \a n can
1543!! be detected by checking .NOT.c_e(this%array(n)).
1544SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1545TYPE(arrayof_georef_coord_array),INTENT(out) :: this !< uninitialised array object
1546CHARACTER(len=*),INTENT(in) :: shpfile !< name of shapefile (with or without extension)
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
1575!> Export an array of \a georef_coord_array objects to a file
1576!! in ESRI/Shapefile format. All the \a this%arraysize shapes
1577!! contained in \a this are exported to the requested
1578!! shapefile. Topology information and possible polygon parts are
1579!! exported as well, while projection information is ignored.
1580SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1581TYPE(arrayof_georef_coord_array),INTENT(in) :: this !< array object to be exported
1582CHARACTER(len=*),INTENT(in) :: shpfile !< name of shapefile (with or without extension)
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
1607!> Determines whether the point \a this lies inside the polygon \a poly.
1608!! The polygon is forced to be closed if it is not already the case,
1609!! and there is no check about the topology of \a poly to really be of
1610!! polygon type. It works also with polygons in parts (as from
1611!! shapefile specification) defining either multiple polygons or
1612!! polygons with holes.
1613!!
1614!! The method used consists in counting the number of intersections as
1615!! indicated in comp.graphics.algorithms FAQ
1616!! (http://www.faqs.org/faqs/graphics/algorithms-faq/) or in
1617!! http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html.
1618FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1619TYPE(georef_coord), INTENT(IN) :: this !< oggetto di cui determinare la posizione
1620TYPE(georef_coord_array), INTENT(IN) :: poly !< poligono contenitore
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 |