libsim Versione 7.1.11

◆ geo_coord_inside()

logical function geo_coord_inside ( type(geo_coord), intent(in)  this,
type(geo_coordvect), intent(in)  poly 
)

Determina se il punto indicato da this si trova dentro o fuori dal poligono descritto dall'oggetto poly.

Funziona anche se la topologia di poly non รจ poligonale, forzandone la chiusura; usa un algoritmo di ricerca del numero di intersezioni, come indicato in comp.graphics.algorithms FAQ (http://www.faqs.org/faqs/graphics/algorithms-faq/) o in http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

Parametri
[in]thisoggetto di cui determinare la posizione
[in]polypoligono contenitore

Definizione alla linea 1055 del file geo_coord_class.F90.

1056! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1057! authors:
1058! Davide Cesari <dcesari@arpa.emr.it>
1059! Paolo Patruno <ppatruno@arpa.emr.it>
1060
1061! This program is free software; you can redistribute it and/or
1062! modify it under the terms of the GNU General Public License as
1063! published by the Free Software Foundation; either version 2 of
1064! the License, or (at your option) any later version.
1065
1066! This program is distributed in the hope that it will be useful,
1067! but WITHOUT ANY WARRANTY; without even the implied warranty of
1068! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1069! GNU General Public License for more details.
1070
1071! You should have received a copy of the GNU General Public License
1072! along with this program. If not, see <http://www.gnu.org/licenses/>.
1073#include "config.h"
1074
1083MODULE geo_coord_class
1084USE kinds
1085USE err_handling
1088!USE doubleprecision_phys_const
1090#ifdef HAVE_SHAPELIB
1091USE shapelib
1092!, ONLY: shpobject, shpreadobject, shpdestroyobject, shpwriteobject, &
1093! shpcreatesimpleobject
1094#endif
1095IMPLICIT NONE
1096
1097
1105INTEGER, PARAMETER :: fp_geo=fp_d
1106real(kind=fp_geo), PARAMETER :: fp_geo_miss=dmiss
1107
1110TYPE geo_coord
1111 PRIVATE
1112 INTEGER(kind=int_l) :: ilon, ilat
1113END TYPE geo_coord
1114
1115TYPE geo_coord_r
1116 PRIVATE
1117 REAL(kind=fp_geo) :: lon, lat
1118END TYPE geo_coord_r
1119
1120
1124TYPE geo_coordvect
1125 PRIVATE
1126 REAL(kind=fp_geo),POINTER :: ll(:,:) => null()
1127 INTEGER :: vsize, vtype
1128END TYPE geo_coordvect
1129
1131TYPE(geo_coord),PARAMETER :: geo_coord_miss=geo_coord(imiss,imiss)
1132
1133INTEGER, PARAMETER :: geo_coordvect_point = 1
1134INTEGER, PARAMETER :: geo_coordvect_arc = 3
1135INTEGER, PARAMETER :: geo_coordvect_polygon = 5
1136INTEGER, PARAMETER :: geo_coordvect_multipoint = 8
1137
1138
1142INTERFACE init
1143 MODULE PROCEDURE geo_coord_init, geo_coordvect_init
1144END INTERFACE
1145
1149INTERFACE delete
1150 MODULE PROCEDURE geo_coord_delete, geo_coordvect_delete
1151END INTERFACE
1152
1154INTERFACE getval
1155 MODULE PROCEDURE geo_coord_getval, geo_coordvect_getval
1156END INTERFACE
1157
1159INTERFACE OPERATOR (==)
1160 MODULE PROCEDURE geo_coord_eq
1161END INTERFACE
1162
1164INTERFACE OPERATOR (/=)
1165 MODULE PROCEDURE geo_coord_ne
1166END INTERFACE
1167
1168INTERFACE OPERATOR (>=)
1169 MODULE PROCEDURE geo_coord_ge
1170END INTERFACE
1171
1172INTERFACE OPERATOR (>)
1173 MODULE PROCEDURE geo_coord_gt
1174END INTERFACE
1175
1176INTERFACE OPERATOR (<=)
1177 MODULE PROCEDURE geo_coord_le
1178END INTERFACE
1179
1180INTERFACE OPERATOR (<)
1181 MODULE PROCEDURE geo_coord_lt
1182END INTERFACE
1183
1186INTERFACE import
1187 MODULE PROCEDURE geo_coordvect_import, geo_coordvect_importvect
1188END INTERFACE
1189
1192INTERFACE export
1193 MODULE PROCEDURE geo_coordvect_export, geo_coordvect_exportvect
1194END INTERFACE
1195
1198INTERFACE read_unit
1199 MODULE PROCEDURE geo_coord_read_unit, geo_coord_vect_read_unit
1200END INTERFACE
1201
1204INTERFACE write_unit
1205 MODULE PROCEDURE geo_coord_write_unit, geo_coord_vect_write_unit
1206END INTERFACE
1207
1209INTERFACE inside
1210 MODULE PROCEDURE geo_coord_inside, geo_coord_inside_rectang
1211END INTERFACE
1212
1214INTERFACE c_e
1215 MODULE PROCEDURE c_e_geo_coord
1216END INTERFACE
1217
1219INTERFACE to_char
1220 MODULE PROCEDURE to_char_geo_coord
1221END INTERFACE
1222
1224INTERFACE display
1225 MODULE PROCEDURE display_geo_coord
1226END INTERFACE
1227
1228CONTAINS
1229
1230
1231! ===================
1232! == geo_coord ==
1233! ===================
1240SUBROUTINE geo_coord_init(this, lon, lat, ilon, ilat)
1241TYPE(geo_coord) :: this
1242REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon
1243REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat
1244integer(kind=int_l), INTENT(IN), OPTIONAL :: ilon
1245integer(kind=int_l), INTENT(IN), OPTIONAL :: ilat
1246
1247real(kind=fp_geo) :: llon,llat
1248
1249CALL optio(ilon, this%ilon)
1250CALL optio(ilat, this%ilat)
1251
1252if (.not. c_e(this%ilon)) then
1253 CALL optio(lon, llon)
1254 if (c_e(llon)) this%ilon=nint(llon*1.d5)
1255end if
1256
1257if (.not. c_e(this%ilat)) then
1258 CALL optio(lat, llat)
1259 if (c_e(llat)) this%ilat=nint(llat*1.d5)
1260end if
1261
1262END SUBROUTINE geo_coord_init
1263
1265SUBROUTINE geo_coord_delete(this)
1266TYPE(geo_coord), INTENT(INOUT) :: this
1267
1268this%ilon = imiss
1269this%ilat = imiss
1270
1271END SUBROUTINE geo_coord_delete
1272
1277elemental SUBROUTINE geo_coord_getval(this, lon, lat,ilon,ilat)
1278TYPE(geo_coord),INTENT(IN) :: this
1279REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lon
1280REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lat
1281integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilon
1282integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilat
1283
1284IF (PRESENT(ilon)) ilon = getilon(this)
1285IF (PRESENT(ilat)) ilat = getilat(this)
1286
1287IF (PRESENT(lon)) lon = getlon(this)
1288IF (PRESENT(lat)) lat = getlat(this)
1289
1290END SUBROUTINE geo_coord_getval
1291
1292
1297elemental FUNCTION getilat(this)
1298TYPE(geo_coord),INTENT(IN) :: this
1299integer(kind=int_l) :: getilat
1300
1301getilat = this%ilat
1302
1303END FUNCTION getilat
1304
1309elemental FUNCTION getlat(this)
1310TYPE(geo_coord),INTENT(IN) :: this
1311real(kind=fp_geo) :: getlat
1312integer(kind=int_l) :: ilat
1313
1314ilat=getilat(this)
1315if (c_e(ilat)) then
1316 getlat = ilat*1.d-5
1317else
1318 getlat=fp_geo_miss
1319end if
1320
1321END FUNCTION getlat
1322
1327elemental FUNCTION getilon(this)
1328TYPE(geo_coord),INTENT(IN) :: this
1329integer(kind=int_l) :: getilon
1330
1331getilon = this%ilon
1332
1333END FUNCTION getilon
1334
1335
1340elemental FUNCTION getlon(this)
1341TYPE(geo_coord),INTENT(IN) :: this
1342real(kind=fp_geo) :: getlon
1343integer(kind=int_l) :: ilon
1344
1345ilon=getilon(this)
1346if (c_e(ilon)) then
1347 getlon = ilon*1.d-5
1348else
1349 getlon=fp_geo_miss
1350end if
1351
1352END FUNCTION getlon
1353
1354
1355elemental FUNCTION geo_coord_eq(this, that) RESULT(res)
1356TYPE(geo_coord),INTENT(IN) :: this, that
1357LOGICAL :: res
1358
1359res = (this%ilon == that%ilon .AND. this%ilat == that%ilat)
1360
1361END FUNCTION geo_coord_eq
1362
1363
1364elemental FUNCTION geo_coord_ne(this, that) RESULT(res)
1365TYPE(geo_coord),INTENT(IN) :: this, that
1366LOGICAL :: res
1367
1368res = .not. this == that
1369
1370END FUNCTION geo_coord_ne
1371
1374elemental FUNCTION geo_coord_gt(this, that) RESULT(res)
1375TYPE(geo_coord),INTENT(IN) :: this, that
1376LOGICAL :: res
1377
1378res = this%ilon > that%ilon
1379
1380if ( this%ilon == that%ilon ) then
1381 res = this%ilat > that%ilat
1382end if
1383
1384END FUNCTION geo_coord_gt
1385
1386
1389elemental FUNCTION geo_coord_ge(this, that) RESULT(res)
1390TYPE(geo_coord),INTENT(IN) :: this, that
1391LOGICAL :: res
1392
1393res = .not. this < that
1394
1395END FUNCTION geo_coord_ge
1396
1399elemental FUNCTION geo_coord_lt(this, that) RESULT(res)
1400TYPE(geo_coord),INTENT(IN) :: this, that
1401LOGICAL :: res
1402
1403res = this%ilon < that%ilon
1404
1405if ( this%ilon == that%ilon ) then
1406 res = this%ilat < that%ilat
1407end if
1408
1409
1410END FUNCTION geo_coord_lt
1411
1412
1415elemental FUNCTION geo_coord_le(this, that) RESULT(res)
1416TYPE(geo_coord),INTENT(IN) :: this, that
1417LOGICAL :: res
1418
1419res = .not. this > that
1420
1421END FUNCTION geo_coord_le
1422
1423
1426elemental FUNCTION geo_coord_ure(this, that) RESULT(res)
1427TYPE(geo_coord),INTENT(IN) :: this, that
1428LOGICAL :: res
1429
1430res = (this%ilon >= that%ilon .AND. this%ilat >= that%ilat)
1431
1432END FUNCTION geo_coord_ure
1433
1436elemental FUNCTION geo_coord_ur(this, that) RESULT(res)
1437TYPE(geo_coord),INTENT(IN) :: this, that
1438LOGICAL :: res
1439
1440res = (this%ilon > that%ilon .AND. this%ilat > that%ilat)
1441
1442END FUNCTION geo_coord_ur
1443
1444
1447elemental FUNCTION geo_coord_lle(this, that) RESULT(res)
1448TYPE(geo_coord),INTENT(IN) :: this, that
1449LOGICAL :: res
1450
1451res = (this%ilon <= that%ilon .AND. this%ilat <= that%ilat)
1452
1453END FUNCTION geo_coord_lle
1454
1457elemental FUNCTION geo_coord_ll(this, that) RESULT(res)
1458TYPE(geo_coord),INTENT(IN) :: this, that
1459LOGICAL :: res
1460
1461res = (this%ilon < that%ilon .AND. this%ilat < that%ilat)
1462
1463END FUNCTION geo_coord_ll
1464
1465
1471SUBROUTINE geo_coord_read_unit(this, unit)
1472TYPE(geo_coord),INTENT(out) :: this
1473INTEGER, INTENT(in) :: unit
1474
1475CALL geo_coord_vect_read_unit((/this/), unit)
1476
1477END SUBROUTINE geo_coord_read_unit
1478
1479
1485SUBROUTINE geo_coord_vect_read_unit(this, unit)
1486TYPE(geo_coord) :: this(:)
1487INTEGER, INTENT(in) :: unit
1488
1489CHARACTER(len=40) :: form
1490INTEGER :: i
1491
1492INQUIRE(unit, form=form)
1493IF (form == 'FORMATTED') THEN
1494 read(unit,*) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1495!TODO bug gfortran compiler !
1496!missing values are unredeable when formatted
1497ELSE
1498 READ(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1499ENDIF
1500
1501END SUBROUTINE geo_coord_vect_read_unit
1502
1503
1508SUBROUTINE geo_coord_write_unit(this, unit)
1509TYPE(geo_coord),INTENT(in) :: this
1510INTEGER, INTENT(in) :: unit
1511
1512CALL geo_coord_vect_write_unit((/this/), unit)
1513
1514END SUBROUTINE geo_coord_write_unit
1515
1516
1521SUBROUTINE geo_coord_vect_write_unit(this, unit)
1522TYPE(geo_coord),INTENT(in) :: this(:)
1523INTEGER, INTENT(in) :: unit
1524
1525CHARACTER(len=40) :: form
1526INTEGER :: i
1527
1528INQUIRE(unit, form=form)
1529IF (form == 'FORMATTED') THEN
1530 WRITE(unit,*) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1531!TODO bug gfortran compiler !
1532!missing values are unredeable when formatted
1533ELSE
1534 WRITE(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1535ENDIF
1536
1537END SUBROUTINE geo_coord_vect_write_unit
1538
1539
1542FUNCTION geo_coord_dist(this, that) RESULT(dist)
1544TYPE(geo_coord), INTENT (IN) :: this
1545TYPE(geo_coord), INTENT (IN) :: that
1546REAL(kind=fp_geo) :: dist
1547
1548REAL(kind=fp_geo) :: x,y
1549! Distanza approssimata, valida per piccoli angoli
1550
1551x = (getlon(this)-getlon(that))*cos(((getlat(this)+getlat(that))/2.)*degrad)
1552y = getlat(this)-getlat(that)
1553dist = sqrt(x**2 + y**2)*degrad*rearth
1554
1555END FUNCTION geo_coord_dist
1556
1557
1566FUNCTION geo_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1567TYPE(geo_coord),INTENT(IN) :: this
1568TYPE(geo_coord),INTENT(IN) :: coordmin
1569TYPE(geo_coord),INTENT(IN) :: coordmax
1570LOGICAL :: res
1571
1572res = (geo_coord_ure(this, coordmin) .AND. geo_coord_lle(this,coordmax))
1573
1574END FUNCTION geo_coord_inside_rectang
1575
1576
1577! ===================
1578! == geo_coordvect ==
1579! ===================
1588RECURSIVE SUBROUTINE geo_coordvect_init(this, lon, lat)
1589TYPE(geo_coordvect), INTENT(OUT) :: this
1590REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon(:)
1591REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat(:)
1592
1593IF (PRESENT(lon) .AND. PRESENT(lat)) THEN
1594 this%vsize = min(SIZE(lon), SIZE(lat))
1595 ALLOCATE(this%ll(this%vsize,2))
1596 this%ll(1:this%vsize,1) = lon(1:this%vsize)
1597 this%ll(1:this%vsize,2) = lat(1:this%vsize)
1598ELSE
1599 this%vsize = 0
1600 NULLIFY(this%ll)
1601ENDIF
1602this%vtype = 0 !?
1603
1604END SUBROUTINE geo_coordvect_init
1605
1606
1609SUBROUTINE geo_coordvect_delete(this)
1610TYPE(geo_coordvect), INTENT(INOUT) :: this
1611
1612IF (ASSOCIATED(this%ll)) DEALLOCATE(this%ll)
1613this%vsize = 0
1614this%vtype = 0
1615
1616END SUBROUTINE geo_coordvect_delete
1617
1618
1627SUBROUTINE geo_coordvect_getval(this, lon, lat)
1628TYPE(geo_coordvect),INTENT(IN) :: this
1629REAL(kind=fp_geo), OPTIONAL, POINTER :: lon(:)
1630REAL(kind=fp_geo), OPTIONAL, POINTER :: lat(:)
1631
1632IF (PRESENT(lon)) THEN
1633 IF (ASSOCIATED(this%ll)) THEN
1634 ALLOCATE(lon(this%vsize))
1635 lon(:) = this%ll(1:this%vsize,1)
1636 ENDIF
1637ENDIF
1638IF (PRESENT(lat)) THEN
1639 IF (ASSOCIATED(this%ll)) THEN
1640 ALLOCATE(lat(this%vsize))
1641 lat(:) = this%ll(1:this%vsize,2)
1642 ENDIF
1643ENDIF
1644
1645END SUBROUTINE geo_coordvect_getval
1646
1647
1648SUBROUTINE geo_coordvect_import(this, unitsim, &
1649#ifdef HAVE_SHAPELIB
1650 shphandle, &
1651#endif
1652 nshp)
1653TYPE(geo_coordvect), INTENT(OUT) :: this
1654INTEGER,OPTIONAL,INTENT(IN) :: unitsim
1655#ifdef HAVE_SHAPELIB
1656TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
1657#endif
1658INTEGER,OPTIONAL,INTENT(IN) :: nshp
1659
1660REAL(kind=fp_geo),ALLOCATABLE :: llon(:), llat(:)
1661REAL(kind=fp_geo) :: lv1,lv2,lv3,lv4
1662INTEGER :: i, lvsize
1663CHARACTER(len=40) :: lname
1664#ifdef HAVE_SHAPELIB
1665TYPE(shpobject) :: shpobj
1666#endif
1667
1668IF (PRESENT(unitsim)) THEN
1669 ! Leggo l'intestazione
1670 READ(unitsim,*,END=10)lvsize,lv1,lv2,lv3,lv4,lname
1671 ALLOCATE(llon(lvsize+1), llat(lvsize+1))
1672 ! Leggo il poligono
1673 READ(unitsim,*)(llon(i),llat(i), i=1,lvsize)
1674 ! Lo chiudo se necessario
1675 IF (llon(1) /= llon(lvsize) .OR. llat(1) /= llat(lvsize)) THEN
1676 lvsize = lvsize + 1
1677 llon(lvsize) = llon(1)
1678 llat(lvsize) = llat(1)
1679 ENDIF
1680 ! Lo inserisco nel mio oggetto
1681 CALL init(this, lon=llon(1:lvsize), lat=llat(1:lvsize))
1682 this%vtype = geo_coordvect_polygon ! Sempre un poligono
1683
1684 DEALLOCATE(llon, llat)
1685 RETURN
168610 CALL raise_error('nella lettura del file '//trim(to_char(unitsim)))
1687 DEALLOCATE(llon, llat) ! End of file, ritorno un oggetto non assegnato
1688#ifdef HAVE_SHAPELIB
1689ELSE IF (PRESENT(shphandle) .AND. PRESENT(nshp)) THEN
1690 ! Leggo l'oggetto shape
1691 shpobj = shpreadobject(shphandle, nshp)
1692 IF (.NOT.shpisnull(shpobj)) THEN
1693 ! Lo inserisco nel mio oggetto
1694 CALL init(this, lon=real(shpobj%padfx,kind=fp_geo), &
1695 lat=real(shpobj%padfy,kind=fp_geo))
1696 this%vtype = shpobj%nshptype
1697 CALL shpdestroyobject(shpobj)
1698 ELSE
1699 CALL init(this)
1700 ENDIF
1701#endif
1702ENDIF
1703
1704END SUBROUTINE geo_coordvect_import
1705
1706
1707SUBROUTINE geo_coordvect_export(this, unitsim, &
1708#ifdef HAVE_SHAPELIB
1709 shphandle, &
1710#endif
1711 nshp)
1712TYPE(geo_coordvect), INTENT(INOUT) :: this
1713INTEGER,OPTIONAL,INTENT(IN) :: unitsim
1714#ifdef HAVE_SHAPELIB
1715TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
1716#endif
1717INTEGER,OPTIONAL,INTENT(IN) :: nshp
1718
1719INTEGER :: i, lnshp
1720#ifdef HAVE_SHAPELIB
1721TYPE(shpobject) :: shpobj
1722#endif
1723
1724IF (PRESENT(unitsim)) THEN
1725 IF (this%vsize > 0) THEN
1726 ! Scrivo l'intestazione
1727 WRITE(unitsim,*)SIZE(this%ll,1),-1.,5000.,-0.1,1.1,'Area'
1728 ! Scrivo il poligono
1729 WRITE(unitsim,*)(this%ll(i,1:2), i=1,this%vsize)
1730 ELSE
1731 CALL raise_warning('oggetto geo_coordvect vuoto, non scrivo niente in '// &
1732 trim(to_char(unitsim)))
1733 ENDIF
1734#ifdef HAVE_SHAPELIB
1735ELSE IF (PRESENT(shphandle)) THEN
1736 IF (PRESENT(nshp)) THEN
1737 lnshp = nshp
1738 ELSE
1739 lnshp = -1 ! -1 = append
1740 ENDIF
1741 ! Creo l'oggetto shape inizializzandolo con il mio oggetto
1742 shpobj = shpcreatesimpleobject(this%vtype, this%vsize, &
1743 REAL(this%ll(1:this%vsize,1),kind=fp_d), &
1744 REAL(this%ll(1:this%vsize,2),kind=fp_d))
1745 IF (.NOT.shpisnull(shpobj)) THEN
1746 ! Lo scrivo nello shapefile
1747 i=shpwriteobject(shphandle, lnshp, shpobj)
1748 CALL shpdestroyobject(shpobj)
1749 ENDIF
1750#endif
1751ENDIF
1752
1753END SUBROUTINE geo_coordvect_export
1754
1773SUBROUTINE geo_coordvect_importvect(this, shpfilesim, shpfile)
1774TYPE(geo_coordvect),POINTER :: this(:)
1775CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim
1776CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile
1777
1778REAL(kind=fp_geo) :: inu
1779REAL(kind=fp_d) :: minb(4), maxb(4)
1780INTEGER :: i, u, ns, lvsize, shptype, dbfnf, dbfnr
1781CHARACTER(len=40) :: lname
1782#ifdef HAVE_SHAPELIB
1783TYPE(shpfileobject) :: shphandle
1784#endif
1785
1786NULLIFY(this)
1787
1788IF (PRESENT(shpfilesim)) THEN
1789 u = getunit()
1790 OPEN(u, file=shpfilesim, status='old', err=30)
1791 ns = 0 ! Conto il numero di shape contenute
1792 DO WHILE(.true.)
1793 READ(u,*,END=10,ERR=20)lvsize,inu,inu,inu,inu,lname
1794 READ(u,*,END=20,ERR=20)(inu,inu, i=1,lvsize)
1795 ns = ns + 1
1796 ENDDO
179710 CONTINUE
1798 IF (ns > 0) THEN ! Alloco e leggo il mio oggetto
1799 ALLOCATE(this(ns))
1800 rewind(u)
1801 DO i = 1, ns
1802 CALL import(this(i), unitsim=u)
1803 ENDDO
1804 ENDIF
180520 CONTINUE
1806 CLOSE(u)
1807 IF (.NOT.ASSOCIATED(this)) THEN
1808 CALL raise_warning('file '//trim(shpfilesim)//' vuoto o corrotto')
1809 ENDIF
1810 RETURN
181130 CONTINUE
1812 CALL raise_error('Impossibile aprire il file '//trim(shpfile))
1813 RETURN
1814
1815ELSE IF (PRESENT(shpfile)) THEN
1816#ifdef HAVE_SHAPELIB
1817 shphandle = shpopen(trim(shpfile), 'rb')
1818 IF (shpfileisnull(shphandle)) THEN
1819 CALL raise_error('Impossibile aprire lo shapefile '//trim(shpfile))
1820 RETURN
1821 ENDIF
1822 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr) ! Ottengo le info sul file
1823 IF (ns > 0) THEN ! Alloco e leggo il mio oggetto
1824 ALLOCATE(this(ns))
1825 this(:)%vtype = shptype
1826 DO i = 1, ns
1827 CALL import(this(i), shphandle=shphandle, nshp=i-1)
1828 ENDDO
1829 ENDIF
1830 CALL shpclose(shphandle)
1831 RETURN
1832#endif
1833ENDIF
1834
1835END SUBROUTINE geo_coordvect_importvect
1836
1837
1840SUBROUTINE geo_coordvect_exportvect(this, shpfilesim, shpfile, append)
1841TYPE(geo_coordvect) :: this(:)
1842CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim
1843CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile
1849LOGICAL,INTENT(in),OPTIONAL :: append
1850
1851REAL(kind=fp_d) :: minb(4), maxb(4)
1852INTEGER :: i, u, ns, shptype, dbfnf, dbfnr
1853LOGICAL :: lappend
1854#ifdef HAVE_SHAPELIB
1855TYPE(shpfileobject) :: shphandle
1856#endif
1857
1858IF (PRESENT(append)) THEN
1859 lappend = append
1860ELSE
1861 lappend = .false.
1862ENDIF
1863IF (PRESENT(shpfilesim)) THEN
1864 u = getunit()
1865 IF (lappend) THEN
1866 OPEN(u, file=shpfilesim, status='unknown', position='append', err=30)
1867 ELSE
1868 OPEN(u, file=shpfilesim, status='unknown', err=30)
1869 ENDIF
1870 DO i = 1, SIZE(this)
1871 CALL export(this(i), unitsim=u)
1872 ENDDO
1873 CLOSE(u)
1874 RETURN
187530 CONTINUE
1876 CALL raise_error('Impossibile aprire il file '//trim(shpfile))
1877 RETURN
1878ELSE IF (PRESENT(shpfile)) THEN
1879#ifdef HAVE_SHAPELIB
1880 IF (lappend) THEN
1881 shphandle = shpopen(trim(shpfile), 'r+b')
1882 IF (shpfileisnull(shphandle)) THEN ! shpopen funziona solo su file esistenti
1883 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
1884 ENDIF
1885 ELSE
1886 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
1887 ENDIF
1888 IF (shpfileisnull(shphandle)) THEN
1889 CALL raise_error('Impossibile aprire lo shapefile '//trim(shpfile))
1890 RETURN
1891 ENDIF
1892 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr) ! Ottengo le info sul file
1893 DO i = 1, SIZE(this)
1894 IF (i > ns .OR. lappend) THEN ! Append shape
1895 CALL export(this(i), shphandle=shphandle)
1896 ELSE ! Overwrite shape
1897 CALL export(this(i), shphandle=shphandle, nshp=i-1)
1898 ENDIF
1899 ENDDO
1900 CALL shpclose(shphandle)
1901 RETURN
1902#endif
1903ENDIF
1904
1905END SUBROUTINE geo_coordvect_exportvect
1906
1907
1916FUNCTION geo_coord_inside(this, poly) RESULT(inside)
1917TYPE(geo_coord), INTENT(IN) :: this
1918TYPE(geo_coordvect), INTENT(IN) :: poly
1919LOGICAL :: inside
1920
1921INTEGER :: i, j, starti
1922
1923inside = .false.
1924IF (all(poly%ll(1,:) == poly%ll(poly%vsize,:))) THEN ! Poligono chiuso
1925 starti = 2
1926 j = 1
1927ELSE ! Poligono non chiuso
1928 starti = 1
1929 j = poly%vsize
1930ENDIF
1931DO i = starti, poly%vsize
1932 IF ((poly%ll(i,2) <= getlat(this) .AND. &
1933 getlat(this) < poly%ll(j,2)) .OR. &
1934 (poly%ll(j,2) <= getlat(this) .AND. &
1935 getlat(this) < poly%ll(i,2))) THEN
1936 IF (getlon(this) < (poly%ll(j,1) - poly%ll(i,1)) * &
1937 (getlat(this) - poly%ll(i,2)) / &
1938 (poly%ll(j,2) - poly%ll(i,2)) + poly%ll(i,1)) THEN
1939 inside = .NOT. inside
1940 ENDIF
1941 ENDIF
1942 j = i
1943ENDDO
1944
1945END FUNCTION geo_coord_inside
1946
1947
1948ELEMENTAL FUNCTION c_e_geo_coord(this) result (res)
1949TYPE(geo_coord),INTENT(in) :: this
1950LOGICAL :: res
1951
1952res = .not. this == geo_coord_miss
1953
1954end FUNCTION c_e_geo_coord
1955
1956
1957character(len=80) function to_char_geo_coord(this)
1958TYPE(geo_coord),INTENT(in) :: this
1959
1960to_char_geo_coord = "GEO_COORD: Lon="// &
1961 trim(to_char(getlon(this),miss="Missing lon",form="(f11.5)"))//&
1962 " Lat="// &
1963 trim(to_char(getlat(this),miss="Missing lat",form="(f11.5)"))
1964
1965end function to_char_geo_coord
1966
1967
1968subroutine display_geo_coord(this)
1969TYPE(geo_coord),INTENT(in) :: this
1970
1971print*,trim(to_char(this))
1972
1973end subroutine display_geo_coord
1974
1975
1976END MODULE geo_coord_class
Detructors for the two classes.
Export one or more geo_coordvect objects to a plain text file or to a file in ESRI/Shapefile format.
Methods for returning the value of object members.
Import one or more geo_coordvect objects from a plain text file or for a file in ESRI/Shapefile forma...
Constructors for the two classes.
Determine whether a point lies inside a polygon or a rectangle.
Read a single geo_coord object or an array of geo_coord objects from a Fortran FORMATTED or UNFORMATT...
Represent geo_coord object in a pretty string.
Write a single geo_coord object or an array of geo_coord objects to a Fortran FORMATTED or UNFORMATTE...
Utilities for CHARACTER variables.
Costanti fisiche (DOUBLEPRECISION).
Definition: phys_const.f90:72
Gestione degli errori.
Utilities for managing files.
Classes for handling georeferenced sparse points in geographical corodinates.
Definition of constants to be used for declaring variables of a desired type.
Definition: kinds.F90:251
Definitions of constants and functions for working with missing values.
Derived type defining an isolated georeferenced point on Earth in polar geographical coordinates.
Derived type defining a one-dimensional array of georeferenced points with an associated topology (is...

Generated with Doxygen.