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