libsim Versione 7.2.6

◆ georef_coord_array_getval()

subroutine georef_coord_array_getval ( type(georef_coord_array), intent(in) this,
double precision, dimension(:), intent(out), optional, allocatable x,
double precision, dimension(:), intent(out), optional, allocatable y,
integer, intent(out), optional topo,
type(geo_proj), intent(out), optional proj )
private

Query a georef_coord_array object.

This is the correct way to retrieve the contents of a georef_coord_array object, since its members are declared as PRIVATE.

Parametri
[in]thisobject to query
[out]xx-coordinate
[out]yy-coordinate
[out]topotopology associated with the coordinates
[out]projgeographical projection

Definizione alla linea 1005 del file georef_coord_class.F90.

1006! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1007! authors:
1008! Davide Cesari <dcesari@arpa.emr.it>
1009! Paolo Patruno <ppatruno@arpa.emr.it>
1010
1011! This program is free software; you can redistribute it and/or
1012! modify it under the terms of the GNU General Public License as
1013! published by the Free Software Foundation; either version 2 of
1014! the License, or (at your option) any later version.
1015
1016! This program is distributed in the hope that it will be useful,
1017! but WITHOUT ANY WARRANTY; without even the implied warranty of
1018! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1019! GNU General Public License for more details.
1020
1021! You should have received a copy of the GNU General Public License
1022! along with this program. If not, see <http://www.gnu.org/licenses/>.
1023#include "config.h"
1024!> This module defines objects describing georeferenced sparse points
1025!! possibly with topology and projection information. This module
1026!! defines two classes, \a georef_coord, which represents a single
1027!! georeferenced point on the Earth, and \a georef_coord_array which
1028!! defines a set of points with a topological relation.
1029!!
1030!! Both classes have \a PRIVATE members, so that they cannot be
1031!! manipulated directly, but only through the proper methods.
1032!!
1033!! It is also possible to dafine a dynamically extendible array of \a
1034!! georef_coord_array objects, of type \a arrayof_georef_coord_array,
1035!! suitable for importing/exporting data from/to a shapefile.
1036!!
1037!! \ingroup base
1038MODULE georef_coord_class
1039USE err_handling
1042USE geo_proj_class
1043#ifdef HAVE_SHAPELIB
1044USE shapelib
1045#endif
1046IMPLICIT NONE
1047
1048!> Derive type defining a single georeferenced point, either in
1049!! geodetic or in projected coordinates. The object has no information
1050!! about the georeferenced coordinate system associated, it must be
1051!! kept separately by the user.
1052TYPE georef_coord
1053 PRIVATE
1054 DOUBLE PRECISION :: x=dmiss, y=dmiss
1055END TYPE georef_coord
1056
1057!> Missing value for georef_coord
1058TYPE(georef_coord),PARAMETER :: georef_coord_miss=georef_coord(dmiss,dmiss)
1059
1060!> Derived type defining a one-dimensional array of georeferenced points
1061!! with an associated topology (isolated point, arc, polygon, group of
1062!! points), possibly broken into parts and with an associated
1063!! georeferenced coordinate system.
1065 PRIVATE
1066 INTEGER,ALLOCATABLE :: parts(:)
1067 TYPE(georef_coord),ALLOCATABLE :: coord(:)
1068 INTEGER :: topo=imiss
1069 TYPE(geo_proj) :: proj
1070 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
1071 LOGICAL :: bbox_updated=.false.
1072END TYPE georef_coord_array
1073
1074INTEGER,PARAMETER :: georef_coord_array_point = 1 !< Topology for georef_coord_array (from shapelib): isolated point
1075INTEGER,PARAMETER :: georef_coord_array_arc = 3 !< Topology for georef_coord_array (from shapelib): arc (multiple arcs unsupported)
1076INTEGER,PARAMETER :: georef_coord_array_polygon = 5 !< Topology for georef_coord_array (from shapelib): polygon (necessarily closed, multiple polygons unsupported)
1077INTEGER,PARAMETER :: georef_coord_array_multipoint = 8 !< Topology for georef_coord_array (from shapelib): group of points
1078
1079
1080!> Detructors for the two classes.
1081!! They clean up all the information associated with the corresponding
1082!! objects.
1083INTERFACE delete
1084 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1085END INTERFACE
1086
1087!> Check missing value.
1088INTERFACE c_e
1089 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1090END INTERFACE
1091
1092!> Methods for returning the value of object members.
1093INTERFACE getval
1094 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
1095END INTERFACE
1096
1097INTERFACE compute_bbox
1098 MODULE PROCEDURE georef_coord_array_compute_bbox
1099END INTERFACE
1100
1101!> Logical equality operator.
1102INTERFACE OPERATOR (==)
1103 MODULE PROCEDURE georef_coord_eq
1104END INTERFACE
1105
1106!> Logical inequality operator.
1107INTERFACE OPERATOR (/=)
1108 MODULE PROCEDURE georef_coord_ne
1109END INTERFACE
1110
1111!> Logical greater-equal operator. Returns true if the first point
1112!! lies to the northeast of the second
1113INTERFACE OPERATOR (>=)
1114 MODULE PROCEDURE georef_coord_ge
1115END INTERFACE
1116
1117!> Logical less-equal operator. Returns true if the first point
1118!! lies to the southwest of the second
1119INTERFACE OPERATOR (<=)
1120 MODULE PROCEDURE georef_coord_le
1121END INTERFACE
1122
1123#ifdef HAVE_SHAPELIB
1124!> Import an array of \a georef_coord_array objects from a file
1125!! in ESRI/Shapefile format.
1126INTERFACE import
1127 MODULE PROCEDURE arrayof_georef_coord_array_import
1128END INTERFACE
1129
1130!> Export an array of \a georef_coord_array objects to a file
1131!! in ESRI/Shapefile format.
1132INTERFACE export
1133 MODULE PROCEDURE arrayof_georef_coord_array_export
1134END INTERFACE
1135#endif
1136
1137!> Read a single \a georef_coord object or an array of \a georef_coord objects
1138!! from a Fortran \c FORMATTED or \c UNFORMATTED file.
1139INTERFACE read_unit
1140 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1141END INTERFACE
1142
1143!> Write a single \a georef_coord object or an array of \a georef_coord objects
1144!! to a Fortran \c FORMATTED or \c UNFORMATTED file.
1145INTERFACE write_unit
1146 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1147END INTERFACE
1148
1149!> Determine whether a point lies inside a polygon or a rectangle.
1150INTERFACE inside
1151 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1152END INTERFACE
1153
1154!> Compute the distance in m between two points
1155INTERFACE dist
1156 MODULE PROCEDURE georef_coord_dist
1157END INTERFACE
1158
1159#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1160#define ARRAYOF_TYPE arrayof_georef_coord_array
1161!define ARRAYOF_ORIGEQ 0
1162#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1163#include "arrayof_pre.F90"
1164! from arrayof
1165PUBLIC insert, append, remove, packarray
1166
1167PRIVATE
1168PUBLIC georef_coord, georef_coord_miss, &
1169 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1170 georef_coord_array_polygon, georef_coord_array_multipoint, &
1171 delete, c_e, getval, compute_bbox, OPERATOR(==), OPERATOR(/=), OPERATOR(>=), OPERATOR(<=), &
1172#ifdef HAVE_SHAPELIB
1173 import, export, &
1174#endif
1176 georef_coord_new, georef_coord_array_new
1177
1178CONTAINS
1179
1180#include "arrayof_post.F90"
1181
1182! ===================
1183! == georef_coord ==
1184! ===================
1185!> Construct a \a georef_coord object with the optional parameters provided.
1186!! If coordinates are not provided the object obtained is empty
1187!! (missing, see c_e function).
1188FUNCTION georef_coord_new(x, y) RESULT(this)
1189DOUBLE PRECISION,INTENT(in),OPTIONAL :: x !< x coordinate
1190DOUBLE PRECISION,INTENT(in),OPTIONAL :: y !< y coordinate
1191TYPE(georef_coord) :: this
1192
1193CALL optio(x, this%x)
1194CALL optio(y, this%y)
1195
1196END FUNCTION georef_coord_new
1197
1198
1199SUBROUTINE georef_coord_delete(this)
1200TYPE(georef_coord),INTENT(inout) :: this
1201
1202this%x = dmiss
1203this%y = dmiss
1204
1205END SUBROUTINE georef_coord_delete
1206
1207
1208ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1209TYPE(georef_coord),INTENT(in) :: this
1210LOGICAL :: res
1211
1212res = .NOT. this == georef_coord_miss
1213
1214END FUNCTION georef_coord_c_e
1215
1216
1217!> Query a \a georef_coord object.
1218!! This is the correct way to retrieve the contents of a \a
1219!! georef_coord object, since its members are declared as \a
1220!! PRIVATE. It is declared as \a ELEMENTAL, thus it works also on
1221!! arrays of any shape and, in that case, the result will hae the same
1222!! shape as \a this.
1223ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1224TYPE(georef_coord),INTENT(in) :: this !< object to query
1225DOUBLE PRECISION,INTENT(out),OPTIONAL :: x !< x-coordinate
1226DOUBLE PRECISION,INTENT(out),OPTIONAL :: y !< y-coordinate
1227
1228IF (PRESENT(x)) x = this%x
1229IF (PRESENT(y)) y = this%y
1230
1231END SUBROUTINE georef_coord_getval
1232
1233
1234!> Query a \a georef_coord object associating a geographical projection to it.
1235!! This method allow to interpret the \a x,y coordinates of a \a
1236!! georef_coord object as projected on a specified geographical
1237!! projection system and retrieve the geodetic longitude and/or
1238!! latitude associated to them. When \a x or \y are requested it works
1239!! as the basic get_val method. It is declared as \a ELEMENTAL, thus
1240!! it works also on arrays of any shape and, in that case, the result
1241!! will hae the same shape as \a this.
1242ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1243TYPE(georef_coord),INTENT(in) :: this !< object to query
1244TYPE(geo_proj),INTENT(in) :: proj !< geographical projection associated to coordinates of \a this
1245DOUBLE PRECISION,INTENT(out),OPTIONAL :: x !< x-coordinate
1246DOUBLE PRECISION,INTENT(out),OPTIONAL :: y !< y-coordinate
1247DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon !< geodetic longitude
1248DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat !< geodetic latitude
1249
1250DOUBLE PRECISION :: llon, llat
1251
1252IF (PRESENT(x)) x = this%x
1253IF (PRESENT(y)) y = this%y
1254IF (PRESENT(lon) .OR. present(lat)) THEN
1255 CALL unproj(proj, this%x, this%y, llon, llat)
1256 IF (PRESENT(lon)) lon = llon
1257 IF (PRESENT(lat)) lat = llat
1258ENDIF
1259
1260END SUBROUTINE georef_coord_proj_getval
1261
1262
1263! document and improve
1264ELEMENTAL FUNCTION getlat(this)
1265TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1266DOUBLE PRECISION :: getlat ! latitudine geografica
1267
1268getlat = this%y ! change!!!
1269
1270END FUNCTION getlat
1271
1272! document and improve
1273ELEMENTAL FUNCTION getlon(this)
1274TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1275DOUBLE PRECISION :: getlon ! longitudine geografica
1276
1277getlon = this%x ! change!!!
1278
1279END FUNCTION getlon
1280
1281
1282ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1283TYPE(georef_coord),INTENT(IN) :: this, that
1284LOGICAL :: res
1285
1286res = (this%x == that%x .AND. this%y == that%y)
1287
1288END FUNCTION georef_coord_eq
1289
1290
1291ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1292TYPE(georef_coord),INTENT(IN) :: this, that
1293LOGICAL :: res
1294
1295res = (this%x >= that%x .AND. this%y >= that%y)
1296
1297END FUNCTION georef_coord_ge
1298
1299
1300ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1301TYPE(georef_coord),INTENT(IN) :: this, that
1302LOGICAL :: res
1303
1304res = (this%x <= that%x .AND. this%y <= that%y)
1305
1306END FUNCTION georef_coord_le
1307
1308
1309ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1310TYPE(georef_coord),INTENT(IN) :: this, that
1311LOGICAL :: res
1312
1313res = .NOT.(this == that)
1314
1315END FUNCTION georef_coord_ne
1316
1317
1318!> Legge da un'unità di file il contenuto dell'oggetto \a this.
1319!! Il record da leggere deve essere stato scritto con la ::write_unit
1320!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
1321!! del vettore devono essere accordate. Il metodo controlla se il file è
1322!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1323SUBROUTINE georef_coord_read_unit(this, unit)
1324TYPE(georef_coord),INTENT(out) :: this !< oggetto da leggere
1325INTEGER, INTENT(in) :: unit !< unità da cui leggere
1326
1327CALL georef_coord_vect_read_unit((/this/), unit)
1328
1329END SUBROUTINE georef_coord_read_unit
1330
1331
1332!> Legge da un'unità di file il contenuto dell'oggetto \a this.
1333!! Il record da leggere deve essere stato scritto con la ::write_unit
1334!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
1335!! del vettore devono essere accordate. Il metodo controlla se il file è
1336!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1337SUBROUTINE georef_coord_vect_read_unit(this, unit)
1338TYPE(georef_coord) :: this(:) !< oggetto da leggere
1339INTEGER, INTENT(in) :: unit !< unità da cui leggere
1340
1341CHARACTER(len=40) :: form
1342INTEGER :: i
1343
1344INQUIRE(unit, form=form)
1345IF (form == 'FORMATTED') THEN
1346 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1347!TODO bug gfortran compiler !
1348!missing values are unredeable when formatted
1349ELSE
1350 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1351ENDIF
1352
1353END SUBROUTINE georef_coord_vect_read_unit
1354
1355
1356!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
1357!! Il record scritto potrà successivamente essere letto con la ::read_unit.
1358!! Il metodo controlla se il file è
1359!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1360SUBROUTINE georef_coord_write_unit(this, unit)
1361TYPE(georef_coord),INTENT(in) :: this !< oggetto da scrivere
1362INTEGER, INTENT(in) :: unit !< unità su cui scrivere
1363
1364CALL georef_coord_vect_write_unit((/this/), unit)
1365
1366END SUBROUTINE georef_coord_write_unit
1367
1368
1369!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
1370!! Il record scritto potrà successivamente essere letto con la ::read_unit.
1371!! Il metodo controlla se il file è
1372!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1373SUBROUTINE georef_coord_vect_write_unit(this, unit)
1374TYPE(georef_coord),INTENT(in) :: this(:) !< oggetto da scrivere
1375INTEGER, INTENT(in) :: unit !< unità su cui scrivere
1376
1377CHARACTER(len=40) :: form
1378INTEGER :: i
1379
1380INQUIRE(unit, form=form)
1381IF (form == 'FORMATTED') THEN
1382 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1383!TODO bug gfortran compiler !
1384!missing values are unredeable when formatted
1385ELSE
1386 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1387ENDIF
1388
1389END SUBROUTINE georef_coord_vect_write_unit
1390
1391
1392!> Restituisce la distanza in m tra 2 oggetti georef_coord.
1393!! La distanza è calcolata approssimativamente ed è valida per piccoli angoli.
1394FUNCTION georef_coord_dist(this, that) RESULT(dist)
1396TYPE(georef_coord), INTENT (IN) :: this !< primo punto
1397TYPE(georef_coord), INTENT (IN) :: that !< secondo punto
1398DOUBLE PRECISION :: dist !< distanza in metri
1399
1400DOUBLE PRECISION :: x,y
1401! Distanza approssimata, valida per piccoli angoli
1402
1403x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1404y = (this%y-that%y)
1405dist = sqrt(x**2 + y**2)*degrad*rearth
1406
1407END FUNCTION georef_coord_dist
1408
1409
1410!> Determines whether the point \a this lies inside a specified rectangle.
1411!! The rectangle is oriented parallely to the coordinate system, its
1412!! lower-left and upper-right vertices are specified by the two
1413!! arguments. The function returns \c .TRUE. also in the case it lies
1414!! exactly on the border of the rectangle.
1415FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1416TYPE(georef_coord),INTENT(IN) :: this !< point to test
1417TYPE(georef_coord),INTENT(IN) :: coordmin !< lower-left vertex
1418TYPE(georef_coord),INTENT(IN) :: coordmax !< upper-right vertex
1419LOGICAL :: res
1420
1421res = (this >= coordmin .AND. this <= coordmax)
1422
1423END FUNCTION georef_coord_inside_rectang
1424
1425
1426! ========================
1427! == georef_coord_array ==
1428! ========================
1429!> Construct a \a georef_coord_array object with the optional parameters provided.
1430!! If coordinates are not provided the object obtained is empty
1431!! (missing, see c_e function), if coordinate arrays are of different
1432!! lengths the \a georef_coord_array is initialised to the shortest
1433!! length provided.
1434FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1435DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:) !< x coordinate array
1436DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:) !< y coordinate array
1437INTEGER,INTENT(in),OPTIONAL :: topo !< topology of the object, one of the \a georef_coord_array_* constants
1438TYPE(geo_proj),INTENT(in),OPTIONAL :: proj !< geographical projection associated to the coordinates
1439TYPE(georef_coord_array) :: this
1440
1441INTEGER :: lsize
1442
1443IF (PRESENT(x) .AND. PRESENT(y)) THEN
1444 lsize = min(SIZE(x), SIZE(y))
1445 ALLOCATE(this%coord(lsize))
1446 this%coord(1:lsize)%x = x(1:lsize)
1447 this%coord(1:lsize)%y = y(1:lsize)
1448ENDIF
1449this%topo = optio_l(topo)
1450IF (PRESENT(proj)) this%proj = proj
1451
1452END FUNCTION georef_coord_array_new
1453
1454
1455SUBROUTINE georef_coord_array_delete(this)
1456TYPE(georef_coord_array),INTENT(inout) :: this
1457
1458TYPE(georef_coord_array) :: lobj
1459
1460this = lobj
1461
1462END SUBROUTINE georef_coord_array_delete
1463
1464
1465ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1466TYPE(georef_coord_array),INTENT(in) :: this
1467LOGICAL :: res
1468
1469res = ALLOCATED(this%coord)
1470
1471END FUNCTION georef_coord_array_c_e
1472
1473
1474!> Query a \a georef_coord_array object.
1475!! This is the correct way to retrieve the contents of a \a
1476!! georef_coord_array object, since its members are declared as \a
1477!! PRIVATE.
1478SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1479TYPE(georef_coord_array),INTENT(in) :: this !< object to query
1480DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:) !< x-coordinate
1481DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:) !< y-coordinate
1482! allocatable per vedere di nascosto l'effetto che fa
1483INTEGER,OPTIONAL,INTENT(out) :: topo !< topology associated with the coordinates
1484TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj !< geographical projection
1485
1486
1487IF (PRESENT(x)) THEN
1488 IF (ALLOCATED(this%coord)) THEN
1489 x = this%coord%x
1490 ENDIF
1491ENDIF
1492IF (PRESENT(y)) THEN
1493 IF (ALLOCATED(this%coord)) THEN
1494 y = this%coord%y
1495 ENDIF
1496ENDIF
1497IF (PRESENT(topo)) topo = this%topo
1498IF (PRESENT(proj)) proj = this%proj ! warning proj has no missing value yet
1499
1500END SUBROUTINE georef_coord_array_getval
1501
1502
1503!> Compute the bounding box of each shape in \a georef_coord_array object.
1504!! The bounding box is computed and stored in the object, it is used
1505!! by the inside() function for speedup; after it is computed the
1506!! object cannot be changed, otherwise the bounding box will not be
1507!! valid.
1508SUBROUTINE georef_coord_array_compute_bbox(this)
1509TYPE(georef_coord_array),INTENT(inout) :: this !< object to manipulate
1510
1511IF (ALLOCATED(this%coord)) THEN
1512 this%bbox(1)%x = minval(this%coord(:)%x)
1513 this%bbox(1)%y = minval(this%coord(:)%y)
1514 this%bbox(2)%x = maxval(this%coord(:)%x)
1515 this%bbox(2)%y = maxval(this%coord(:)%y)
1516 this%bbox_updated = .true.
1517ENDIF
1518
1519END SUBROUTINE georef_coord_array_compute_bbox
1520
1521#ifdef HAVE_SHAPELIB
1522! internal method for importing a single shape
1523SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1524TYPE(georef_coord_array),INTENT(OUT) :: this
1525TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1526INTEGER,INTENT(IN) :: nshp
1527
1528TYPE(shpobject) :: shpobj
1529
1530! read shape object
1531shpobj = shpreadobject(shphandle, nshp)
1532IF (.NOT.shpisnull(shpobj)) THEN
1533! import it in georef_coord object
1534 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1535 topo=shpobj%nshptype)
1536 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1537 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1538 ELSE IF (ALLOCATED(this%parts)) THEN
1539 DEALLOCATE(this%parts)
1540 ENDIF
1541 CALL shpdestroyobject(shpobj)
1542 CALL compute_bbox(this)
1543ENDIF
1544
1545
1546END SUBROUTINE georef_coord_array_import
1547
1548
1549! internal method for exporting a single shape
1550SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1551TYPE(georef_coord_array),INTENT(in) :: this
1552TYPE(shpfileobject),INTENT(inout) :: shphandle
1553INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1554
1555INTEGER :: i
1556TYPE(shpobject) :: shpobj
1557
1558IF (ALLOCATED(this%coord)) THEN
1559 IF (ALLOCATED(this%parts)) THEN
1560 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1561 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1562 ELSE
1563 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1564 this%coord(:)%x, this%coord(:)%y)
1565 ENDIF
1566ELSE
1567 RETURN
1568ENDIF
1569
1570IF (.NOT.shpisnull(shpobj)) THEN
1571 i = shpwriteobject(shphandle, nshp, shpobj)
1572 CALL shpdestroyobject(shpobj)
1573ENDIF
1574
1575END SUBROUTINE georef_coord_array_export
1576
1577
1578!> Import an array of \a georef_coord_array objects from a file
1579!! in ESRI/Shapefile format. The \a this argument is an uninitialised
1580!! \a arrayof_georef_coord_array, every element of which,
1581!! this%array(n), is of type \a georef_coord_array and, on return,
1582!! will contain information from the n-th shape of the file. Topology
1583!! information and possible polygon parts are imported as well, while
1584!! no projection information, even if available, is imported. An
1585!! error condition while opening the file can be detected by checking
1586!! .NOT.ASSOCIATED(this%array), while an error reading shape \a n can
1587!! be detected by checking .NOT.c_e(this%array(n)).
1588SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1589TYPE(arrayof_georef_coord_array),INTENT(out) :: this !< uninitialised array object
1590CHARACTER(len=*),INTENT(in) :: shpfile !< name of shapefile (with or without extension)
1591
1592REAL(kind=fp_d) :: minb(4), maxb(4)
1593INTEGER :: i, ns, shptype, dbfnf, dbfnr
1594TYPE(shpfileobject) :: shphandle
1595
1596shphandle = shpopen(trim(shpfile), 'rb')
1597IF (shpfileisnull(shphandle)) THEN
1598 ! log here
1599 CALL raise_error()
1600 RETURN
1601ENDIF
1602
1603! get info about file
1604CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1605IF (ns > 0) THEN ! allocate and read the object
1606 CALL insert(this, nelem=ns)
1607 DO i = 1, ns
1608 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1609 ENDDO
1610ENDIF
1611
1612CALL shpclose(shphandle)
1613! pack object to save memory
1614CALL packarray(this)
1615
1616END SUBROUTINE arrayof_georef_coord_array_import
1617
1618
1619!> Export an array of \a georef_coord_array objects to a file
1620!! in ESRI/Shapefile format. All the \a this%arraysize shapes
1621!! contained in \a this are exported to the requested
1622!! shapefile. Topology information and possible polygon parts are
1623!! exported as well, while projection information is ignored.
1624SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1625TYPE(arrayof_georef_coord_array),INTENT(in) :: this !< array object to be exported
1626CHARACTER(len=*),INTENT(in) :: shpfile !< name of shapefile (with or without extension)
1627
1628INTEGER :: i
1629TYPE(shpfileobject) :: shphandle
1630
1631IF (this%arraysize > 0) THEN
1632 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1633ELSE
1634 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1635ENDIF
1636IF (shpfileisnull(shphandle)) THEN
1637 ! log here
1638 CALL raise_error()
1639 RETURN
1640ENDIF
1641
1642DO i = 1, this%arraysize
1643 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1644ENDDO
1645
1646CALL shpclose(shphandle)
1647
1648END SUBROUTINE arrayof_georef_coord_array_export
1649#endif
1650
1651!> Determines whether the point \a this lies inside the polygon \a poly.
1652!! The polygon is forced to be closed if it is not already the case,
1653!! and there is no check about the topology of \a poly to really be of
1654!! polygon type. It works also with polygons in parts (as from
1655!! shapefile specification) defining either multiple polygons or
1656!! polygons with holes.
1657!!
1658!! The method used consists in counting the number of intersections as
1659!! indicated in comp.graphics.algorithms FAQ
1660!! (http://www.faqs.org/faqs/graphics/algorithms-faq/) or in
1661!! http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html.
1662FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1663TYPE(georef_coord), INTENT(IN) :: this !< oggetto di cui determinare la posizione
1664TYPE(georef_coord_array), INTENT(IN) :: poly !< poligono contenitore
1665LOGICAL :: inside
1666
1667INTEGER :: i
1668
1669inside = .false.
1670IF (.NOT.c_e(this)) RETURN
1671IF (.NOT.ALLOCATED(poly%coord)) RETURN
1672! if outside bounding box stop here
1673IF (poly%bbox_updated) THEN
1674 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1675ENDIF
1676
1677IF (ALLOCATED(poly%parts)) THEN
1678 DO i = 1, SIZE(poly%parts)-1
1679 inside = inside .NEQV. pointinpoly(this%x, this%y, &
1680 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1681 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1682 ENDDO
1683 IF (SIZE(poly%parts) > 0) THEN ! safety check
1684 inside = inside .NEQV. pointinpoly(this%x, this%y, &
1685 poly%coord(poly%parts(i)+1:)%x, &
1686 poly%coord(poly%parts(i)+1:)%y)
1687 ENDIF
1688
1689ELSE
1690 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1691 inside = pointinpoly(this%x, this%y, &
1692 poly%coord(:)%x, poly%coord(:)%y)
1693ENDIF
1694
1695CONTAINS
1696
1697FUNCTION pointinpoly(x, y, px, py)
1698DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1699LOGICAL :: pointinpoly
1700
1701INTEGER :: i, j, starti
1702
1703pointinpoly = .false.
1704
1705IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1706 starti = 2
1707 j = 1
1708ELSE ! unclosed polygon
1709 starti = 1
1710 j = SIZE(px)
1711ENDIF
1712DO i = starti, SIZE(px)
1713 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1714 (py(j) <= y .AND. y < py(i))) THEN
1715 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1716 pointinpoly = .NOT. pointinpoly
1717 ENDIF
1718 ENDIF
1719 j = i
1720ENDDO
1721
1722END FUNCTION pointinpoly
1723
1724END FUNCTION georef_coord_inside
1725
1726
1727
1728END MODULE georef_coord_class
Compute forward coordinate transformation from geographical system to projected system.
Compute backward coordinate transformation from projected system to geographical system.
Quick method to append an element to the array.
Detructors for the two classes.
Compute the distance in m between two points.
Export an array of georef_coord_array objects to a file in ESRI/Shapefile format.
Methods for returning the value of object members.
Import an array of georef_coord_array objects from a file in ESRI/Shapefile format.
Method for inserting elements of the array at a desired position.
Determine whether a point lies inside a polygon or a rectangle.
Method for packing the array object reducing at a minimum the memory occupation, without destroying i...
Read a single georef_coord object or an array of georef_coord objects from a Fortran FORMATTED or UNF...
Method for removing elements of the array at a desired position.
Write a single georef_coord object or an array of georef_coord objects to a Fortran FORMATTED or UNFO...
Generic subroutine for checking OPTIONAL parameters.
Costanti fisiche (DOUBLEPRECISION).
Gestione degli errori.
This module defines objects describing georeferenced sparse points possibly with topology and project...
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Derived type defining a one-dimensional array of georeferenced points with an associated topology (is...
Derive type defining a single georeferenced point, either in geodetic or in projected coordinates.

Generated with Doxygen.