libsim Versione 7.2.6

◆ georef_coord_inside_rectang()

logical function georef_coord_inside_rectang ( type(georef_coord), intent(in) this,
type(georef_coord), intent(in) coordmin,
type(georef_coord), intent(in) coordmax )

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 .TRUE. also in the case it lies exactly on the border of the rectangle.

Parametri
[in]thispoint to test
[in]coordminlower-left vertex
[in]coordmaxupper-right vertex

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!> This module defines objects describing georeferenced sparse points
962!! possibly with topology and projection information. This module
963!! defines two classes, \a georef_coord, which represents a single
964!! georeferenced point on the Earth, and \a georef_coord_array which
965!! defines a set of points with a topological relation.
966!!
967!! Both classes have \a PRIVATE members, so that they cannot be
968!! manipulated directly, but only through the proper methods.
969!!
970!! It is also possible to dafine a dynamically extendible array of \a
971!! georef_coord_array objects, of type \a arrayof_georef_coord_array,
972!! suitable for importing/exporting data from/to a shapefile.
973!!
974!! \ingroup base
976USE err_handling
979USE geo_proj_class
980#ifdef HAVE_SHAPELIB
981USE shapelib
982#endif
983IMPLICIT NONE
984
985!> Derive type defining a single georeferenced point, either in
986!! geodetic or in projected coordinates. The object has no information
987!! about the georeferenced coordinate system associated, it must be
988!! kept separately by the user.
989TYPE georef_coord
990 PRIVATE
991 DOUBLE PRECISION :: x=dmiss, y=dmiss
992END TYPE georef_coord
993
994!> Missing value for georef_coord
995TYPE(georef_coord),PARAMETER :: georef_coord_miss=georef_coord(dmiss,dmiss)
996
997!> Derived type defining a one-dimensional array of georeferenced points
998!! with an associated topology (isolated point, arc, polygon, group of
999!! points), possibly broken into parts and with an associated
1000!! georeferenced coordinate system.
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.
1009END TYPE georef_coord_array
1010
1011INTEGER,PARAMETER :: georef_coord_array_point = 1 !< Topology for georef_coord_array (from shapelib): isolated point
1012INTEGER,PARAMETER :: georef_coord_array_arc = 3 !< Topology for georef_coord_array (from shapelib): arc (multiple arcs unsupported)
1013INTEGER,PARAMETER :: georef_coord_array_polygon = 5 !< Topology for georef_coord_array (from shapelib): polygon (necessarily closed, multiple polygons unsupported)
1014INTEGER,PARAMETER :: georef_coord_array_multipoint = 8 !< Topology for georef_coord_array (from shapelib): group of points
1015
1016
1017!> Detructors for the two classes.
1018!! They clean up all the information associated with the corresponding
1019!! objects.
1020INTERFACE delete
1021 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1022END INTERFACE
1023
1024!> Check missing value.
1025INTERFACE c_e
1026 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1027END INTERFACE
1028
1029!> Methods for returning the value of object members.
1030INTERFACE getval
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
1038!> Logical equality operator.
1039INTERFACE OPERATOR (==)
1040 MODULE PROCEDURE georef_coord_eq
1041END INTERFACE
1042
1043!> Logical inequality operator.
1044INTERFACE OPERATOR (/=)
1045 MODULE PROCEDURE georef_coord_ne
1046END INTERFACE
1047
1048!> Logical greater-equal operator. Returns true if the first point
1049!! lies to the northeast of the second
1050INTERFACE OPERATOR (>=)
1051 MODULE PROCEDURE georef_coord_ge
1052END INTERFACE
1053
1054!> Logical less-equal operator. Returns true if the first point
1055!! lies to the southwest of the second
1056INTERFACE OPERATOR (<=)
1057 MODULE PROCEDURE georef_coord_le
1058END INTERFACE
1059
1060#ifdef HAVE_SHAPELIB
1061!> Import an array of \a georef_coord_array objects from a file
1062!! in ESRI/Shapefile format.
1063INTERFACE import
1064 MODULE PROCEDURE arrayof_georef_coord_array_import
1065END INTERFACE
1066
1067!> Export an array of \a georef_coord_array objects to a file
1068!! in ESRI/Shapefile format.
1069INTERFACE export
1070 MODULE PROCEDURE arrayof_georef_coord_array_export
1071END INTERFACE
1072#endif
1073
1074!> Read a single \a georef_coord object or an array of \a georef_coord objects
1075!! from a Fortran \c FORMATTED or \c UNFORMATTED file.
1076INTERFACE read_unit
1077 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1078END INTERFACE
1079
1080!> Write a single \a georef_coord object or an array of \a georef_coord objects
1081!! to a Fortran \c FORMATTED or \c UNFORMATTED file.
1082INTERFACE write_unit
1083 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1084END INTERFACE
1085
1086!> Determine whether a point lies inside a polygon or a rectangle.
1087INTERFACE inside
1088 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1089END INTERFACE
1090
1091!> Compute the distance in m between two points
1092INTERFACE dist
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
1102PUBLIC insert, append, remove, packarray
1103
1104PRIVATE
1105PUBLIC georef_coord, georef_coord_miss, &
1106 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1107 georef_coord_array_polygon, georef_coord_array_multipoint, &
1108 delete, c_e, getval, compute_bbox, OPERATOR(==), OPERATOR(/=), OPERATOR(>=), OPERATOR(<=), &
1109#ifdef HAVE_SHAPELIB
1110 import, export, &
1111#endif
1113 georef_coord_new, georef_coord_array_new
1114
1115CONTAINS
1116
1117#include "arrayof_post.F90"
1118
1119! ===================
1120! == georef_coord ==
1121! ===================
1122!> Construct a \a georef_coord object with the optional parameters provided.
1123!! If coordinates are not provided the object obtained is empty
1124!! (missing, see c_e function).
1125FUNCTION georef_coord_new(x, y) RESULT(this)
1126DOUBLE PRECISION,INTENT(in),OPTIONAL :: x !< x coordinate
1127DOUBLE PRECISION,INTENT(in),OPTIONAL :: y !< y coordinate
1128TYPE(georef_coord) :: this
1129
1130CALL optio(x, this%x)
1131CALL optio(y, this%y)
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
1154!> Query a \a georef_coord object.
1155!! This is the correct way to retrieve the contents of a \a
1156!! georef_coord object, since its members are declared as \a
1157!! PRIVATE. It is declared as \a ELEMENTAL, thus it works also on
1158!! arrays of any shape and, in that case, the result will hae the same
1159!! shape as \a this.
1160ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1161TYPE(georef_coord),INTENT(in) :: this !< object to query
1162DOUBLE PRECISION,INTENT(out),OPTIONAL :: x !< x-coordinate
1163DOUBLE PRECISION,INTENT(out),OPTIONAL :: y !< y-coordinate
1164
1165IF (PRESENT(x)) x = this%x
1166IF (PRESENT(y)) y = this%y
1167
1168END SUBROUTINE georef_coord_getval
1169
1170
1171!> Query a \a georef_coord object associating a geographical projection to it.
1172!! This method allow to interpret the \a x,y coordinates of a \a
1173!! georef_coord object as projected on a specified geographical
1174!! projection system and retrieve the geodetic longitude and/or
1175!! latitude associated to them. When \a x or \y are requested it works
1176!! as the basic get_val method. It is declared as \a ELEMENTAL, thus
1177!! it works also on arrays of any shape and, in that case, the result
1178!! will hae the same shape as \a this.
1179ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1180TYPE(georef_coord),INTENT(in) :: this !< object to query
1181TYPE(geo_proj),INTENT(in) :: proj !< geographical projection associated to coordinates of \a this
1182DOUBLE PRECISION,INTENT(out),OPTIONAL :: x !< x-coordinate
1183DOUBLE PRECISION,INTENT(out),OPTIONAL :: y !< y-coordinate
1184DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon !< geodetic longitude
1185DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat !< geodetic latitude
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
1192 CALL unproj(proj, this%x, this%y, llon, llat)
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
1255!> Legge da un'unità di file il contenuto dell'oggetto \a this.
1256!! Il record da leggere deve essere stato scritto con la ::write_unit
1257!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
1258!! del vettore devono essere accordate. Il metodo controlla se il file è
1259!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1260SUBROUTINE georef_coord_read_unit(this, unit)
1261TYPE(georef_coord),INTENT(out) :: this !< oggetto da leggere
1262INTEGER, INTENT(in) :: unit !< unità da cui leggere
1263
1264CALL georef_coord_vect_read_unit((/this/), unit)
1265
1266END SUBROUTINE georef_coord_read_unit
1267
1268
1269!> Legge da un'unità di file il contenuto dell'oggetto \a this.
1270!! Il record da leggere deve essere stato scritto con la ::write_unit
1271!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
1272!! del vettore devono essere accordate. Il metodo controlla se il file è
1273!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1274SUBROUTINE georef_coord_vect_read_unit(this, unit)
1275TYPE(georef_coord) :: this(:) !< oggetto da leggere
1276INTEGER, INTENT(in) :: unit !< unità da cui leggere
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
1293!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
1294!! Il record scritto potrà successivamente essere letto con la ::read_unit.
1295!! Il metodo controlla se il file è
1296!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1297SUBROUTINE georef_coord_write_unit(this, unit)
1298TYPE(georef_coord),INTENT(in) :: this !< oggetto da scrivere
1299INTEGER, INTENT(in) :: unit !< unità su cui scrivere
1300
1301CALL georef_coord_vect_write_unit((/this/), unit)
1302
1303END SUBROUTINE georef_coord_write_unit
1304
1305
1306!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
1307!! Il record scritto potrà successivamente essere letto con la ::read_unit.
1308!! Il metodo controlla se il file è
1309!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1310SUBROUTINE georef_coord_vect_write_unit(this, unit)
1311TYPE(georef_coord),INTENT(in) :: this(:) !< oggetto da scrivere
1312INTEGER, INTENT(in) :: unit !< unità su cui scrivere
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
1329!> Restituisce la distanza in m tra 2 oggetti georef_coord.
1330!! La distanza è calcolata approssimativamente ed è valida per piccoli angoli.
1331FUNCTION georef_coord_dist(this, that) RESULT(dist)
1333TYPE(georef_coord), INTENT (IN) :: this !< primo punto
1334TYPE(georef_coord), INTENT (IN) :: that !< secondo punto
1335DOUBLE PRECISION :: dist !< distanza in metri
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
1347!> Determines whether the point \a this lies inside a specified rectangle.
1348!! The rectangle is oriented parallely to the coordinate system, its
1349!! lower-left and upper-right vertices are specified by the two
1350!! arguments. The function returns \c .TRUE. also in the case it lies
1351!! exactly on the border of the rectangle.
1352FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1353TYPE(georef_coord),INTENT(IN) :: this !< point to test
1354TYPE(georef_coord),INTENT(IN) :: coordmin !< lower-left vertex
1355TYPE(georef_coord),INTENT(IN) :: coordmax !< upper-right vertex
1356LOGICAL :: res
1357
1358res = (this >= coordmin .AND. this <= coordmax)
1359
1360END FUNCTION georef_coord_inside_rectang
1361
1362
1363! ========================
1364! == georef_coord_array ==
1365! ========================
1366!> Construct a \a georef_coord_array object with the optional parameters provided.
1367!! If coordinates are not provided the object obtained is empty
1368!! (missing, see c_e function), if coordinate arrays are of different
1369!! lengths the \a georef_coord_array is initialised to the shortest
1370!! length provided.
1371FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1372DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:) !< x coordinate array
1373DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:) !< y coordinate array
1374INTEGER,INTENT(in),OPTIONAL :: topo !< topology of the object, one of the \a georef_coord_array_* constants
1375TYPE(geo_proj),INTENT(in),OPTIONAL :: proj !< geographical projection associated to the coordinates
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)
1387IF (PRESENT(proj)) this%proj = proj
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
1411!> Query a \a georef_coord_array object.
1412!! This is the correct way to retrieve the contents of a \a
1413!! georef_coord_array object, since its members are declared as \a
1414!! PRIVATE.
1415SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1416TYPE(georef_coord_array),INTENT(in) :: this !< object to query
1417DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:) !< x-coordinate
1418DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:) !< y-coordinate
1419! allocatable per vedere di nascosto l'effetto che fa
1420INTEGER,OPTIONAL,INTENT(out) :: topo !< topology associated with the coordinates
1421TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj !< geographical projection
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
1435IF (PRESENT(proj)) proj = this%proj ! warning proj has no missing value yet
1436
1437END SUBROUTINE georef_coord_array_getval
1438
1439
1440!> Compute the bounding box of each shape in \a georef_coord_array object.
1441!! The bounding box is computed and stored in the object, it is used
1442!! by the inside() function for speedup; after it is computed the
1443!! object cannot be changed, otherwise the bounding box will not be
1444!! valid.
1445SUBROUTINE georef_coord_array_compute_bbox(this)
1446TYPE(georef_coord_array),INTENT(inout) :: this !< object to manipulate
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
1515!> Import an array of \a georef_coord_array objects from a file
1516!! in ESRI/Shapefile format. The \a this argument is an uninitialised
1517!! \a arrayof_georef_coord_array, every element of which,
1518!! this%array(n), is of type \a georef_coord_array and, on return,
1519!! will contain information from the n-th shape of the file. Topology
1520!! information and possible polygon parts are imported as well, while
1521!! no projection information, even if available, is imported. An
1522!! error condition while opening the file can be detected by checking
1523!! .NOT.ASSOCIATED(this%array), while an error reading shape \a n can
1524!! be detected by checking .NOT.c_e(this%array(n)).
1525SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1526TYPE(arrayof_georef_coord_array),INTENT(out) :: this !< uninitialised array object
1527CHARACTER(len=*),INTENT(in) :: shpfile !< name of shapefile (with or without extension)
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
1543 CALL insert(this, nelem=ns)
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
1551CALL packarray(this)
1552
1553END SUBROUTINE arrayof_georef_coord_array_import
1554
1555
1556!> Export an array of \a georef_coord_array objects to a file
1557!! in ESRI/Shapefile format. All the \a this%arraysize shapes
1558!! contained in \a this are exported to the requested
1559!! shapefile. Topology information and possible polygon parts are
1560!! exported as well, while projection information is ignored.
1561SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1562TYPE(arrayof_georef_coord_array),INTENT(in) :: this !< array object to be exported
1563CHARACTER(len=*),INTENT(in) :: shpfile !< name of shapefile (with or without extension)
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
1588!> Determines whether the point \a this lies inside the polygon \a poly.
1589!! The polygon is forced to be closed if it is not already the case,
1590!! and there is no check about the topology of \a poly to really be of
1591!! polygon type. It works also with polygons in parts (as from
1592!! shapefile specification) defining either multiple polygons or
1593!! polygons with holes.
1594!!
1595!! The method used consists in counting the number of intersections as
1596!! indicated in comp.graphics.algorithms FAQ
1597!! (http://www.faqs.org/faqs/graphics/algorithms-faq/) or in
1598!! http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html.
1599FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1600TYPE(georef_coord), INTENT(IN) :: this !< oggetto di cui determinare la posizione
1601TYPE(georef_coord_array), INTENT(IN) :: poly !< poligono contenitore
1602LOGICAL :: inside
1603
1604INTEGER :: i
1605
1606inside = .false.
1607IF (.NOT.c_e(this)) RETURN
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
1616 inside = inside .NEQV. pointinpoly(this%x, this%y, &
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
1621 inside = inside .NEQV. pointinpoly(this%x, this%y, &
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
1665END 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.