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