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