libsim Versione 7.2.4
|
◆ arrayof_georef_coord_array_import()
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)).
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
1152USE geo_proj_class
1153#ifdef HAVE_SHAPELIB
1154USE shapelib
1155#endif
1156IMPLICIT NONE
1157
1163 PRIVATE
1164 DOUBLE PRECISION :: x=dmiss, y=dmiss
1166
1169
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.
1183
1184INTEGER,PARAMETER :: georef_coord_array_point = 1
1185INTEGER,PARAMETER :: georef_coord_array_arc = 3
1186INTEGER,PARAMETER :: georef_coord_array_polygon = 5
1187INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
1188
1189
1194 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1195END INTERFACE
1196
1199 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1200END INTERFACE
1201
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
1212INTERFACE OPERATOR (==)
1213 MODULE PROCEDURE georef_coord_eq
1214END INTERFACE
1215
1217INTERFACE OPERATOR (/=)
1218 MODULE PROCEDURE georef_coord_ne
1219END INTERFACE
1220
1223INTERFACE OPERATOR (>=)
1224 MODULE PROCEDURE georef_coord_ge
1225END INTERFACE
1226
1229INTERFACE OPERATOR (<=)
1230 MODULE PROCEDURE georef_coord_le
1231END INTERFACE
1232
1233#ifdef HAVE_SHAPELIB
1234
1236INTERFACE import
1237 MODULE PROCEDURE arrayof_georef_coord_array_import
1238END INTERFACE
1239
1243 MODULE PROCEDURE arrayof_georef_coord_array_export
1244END INTERFACE
1245#endif
1246
1250 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1251END INTERFACE
1252
1256 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1257END INTERFACE
1258
1261 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1262END INTERFACE
1263
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
1276
1277PRIVATE
1279 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1280 georef_coord_array_polygon, georef_coord_array_multipoint, &
1282#ifdef HAVE_SHAPELIB
1284#endif
1286 georef_coord_new, georef_coord_array_new
1287
1288CONTAINS
1289
1290#include "arrayof_post.F90"
1291
1292! ===================
1293! == georef_coord ==
1294! ===================
1298FUNCTION georef_coord_new(x, y) RESULT(this)
1299DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1300DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1301TYPE(georef_coord) :: this
1302
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
1333ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1334TYPE(georef_coord),INTENT(in) :: this
1335DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1336DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1337
1338IF (PRESENT(x)) x = this%x
1339IF (PRESENT(y)) y = this%y
1340
1341END SUBROUTINE georef_coord_getval
1342
1343
1352ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1353TYPE(georef_coord),INTENT(in) :: this
1354TYPE(geo_proj),INTENT(in) :: proj
1355DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1356DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1357DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1358DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
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
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
1433SUBROUTINE georef_coord_read_unit(this, unit)
1434TYPE(georef_coord),INTENT(out) :: this
1435INTEGER, INTENT(in) :: unit
1436
1437CALL georef_coord_vect_read_unit((/this/), unit)
1438
1439END SUBROUTINE georef_coord_read_unit
1440
1441
1447SUBROUTINE georef_coord_vect_read_unit(this, unit)
1448TYPE(georef_coord) :: this(:)
1449INTEGER, INTENT(in) :: unit
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
1470SUBROUTINE georef_coord_write_unit(this, unit)
1471TYPE(georef_coord),INTENT(in) :: this
1472INTEGER, INTENT(in) :: unit
1473
1474CALL georef_coord_vect_write_unit((/this/), unit)
1475
1476END SUBROUTINE georef_coord_write_unit
1477
1478
1483SUBROUTINE georef_coord_vect_write_unit(this, unit)
1484TYPE(georef_coord),INTENT(in) :: this(:)
1485INTEGER, INTENT(in) :: unit
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
1504FUNCTION georef_coord_dist(this, that) RESULT(dist)
1506TYPE(georef_coord), INTENT (IN) :: this
1507TYPE(georef_coord), INTENT (IN) :: that
1508DOUBLE PRECISION :: dist
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
1525FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1526TYPE(georef_coord),INTENT(IN) :: this
1527TYPE(georef_coord),INTENT(IN) :: coordmin
1528TYPE(georef_coord),INTENT(IN) :: coordmax
1529LOGICAL :: res
1530
1531res = (this >= coordmin .AND. this <= coordmax)
1532
1533END FUNCTION georef_coord_inside_rectang
1534
1535
1536! ========================
1537! == georef_coord_array ==
1538! ========================
1544FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1545DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1546DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1547INTEGER,INTENT(in),OPTIONAL :: topo
1548TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
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)
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
1588SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1589TYPE(georef_coord_array),INTENT(in) :: this
1590DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1591DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1592! allocatable per vedere di nascosto l'effetto che fa
1593INTEGER,OPTIONAL,INTENT(out) :: topo
1594TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
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
1609
1610END SUBROUTINE georef_coord_array_getval
1611
1612
1618SUBROUTINE georef_coord_array_compute_bbox(this)
1619TYPE(georef_coord_array),INTENT(inout) :: this
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
1698SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1699TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1700CHARACTER(len=*),INTENT(in) :: shpfile
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
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
1725
1726END SUBROUTINE arrayof_georef_coord_array_import
1727
1728
1734SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1735TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1736CHARACTER(len=*),INTENT(in) :: shpfile
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
1772FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1773TYPE(georef_coord), INTENT(IN) :: this
1774TYPE(georef_coord_array), INTENT(IN) :: poly
1775LOGICAL :: inside
1776
1777INTEGER :: i
1778
1779inside = .false.
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
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
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
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 |