libsim Versione 7.2.6

◆ arrayof_georef_coord_array_import()

subroutine arrayof_georef_coord_array_import ( type(arrayof_georef_coord_array), intent(out) this,
character(len=*), intent(in) shpfile )
private

Import an array of georef_coord_array objects from a file in ESRI/Shapefile format.

The this argument is an uninitialised arrayof_georef_coord_array, every element of which, thisarray(n), is of type georef_coord_array and, on return, will contain information from the n-th shape of the file. Topology information and possible polygon parts are imported as well, while no projection information, even if available, is imported. An error condition while opening the file can be detected by checking .NOT.ASSOCIATED(thisarray), while an error reading shape n can be detected by checking .NOT.c_e(thisarray(n)).

Parametri
[out]thisuninitialised array object
[in]shpfilename of shapefile (with or without extension)

Definizione alla linea 1115 del file georef_coord_class.F90.

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