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