libsim Versione 7.2.0

◆ geo_coord_inside()

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

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 1049 del file geo_coord_class.F90.

1050! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1051! authors:
1052! Davide Cesari <dcesari@arpa.emr.it>
1053! Paolo Patruno <ppatruno@arpa.emr.it>
1054
1055! This program is free software; you can redistribute it and/or
1056! modify it under the terms of the GNU General Public License as
1057! published by the Free Software Foundation; either version 2 of
1058! the License, or (at your option) any later version.
1059
1060! This program is distributed in the hope that it will be useful,
1061! but WITHOUT ANY WARRANTY; without even the implied warranty of
1062! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1063! GNU General Public License for more details.
1064
1065! You should have received a copy of the GNU General Public License
1066! along with this program. If not, see <http://www.gnu.org/licenses/>.
1067#include "config.h"
1068
1077MODULE geo_coord_class
1078USE kinds
1079USE err_handling
1082!USE doubleprecision_phys_const
1084#ifdef HAVE_SHAPELIB
1085USE shapelib
1086!, ONLY: shpobject, shpreadobject, shpdestroyobject, shpwriteobject, &
1087! shpcreatesimpleobject
1088#endif
1089IMPLICIT NONE
1090
1091
1099INTEGER, PARAMETER :: fp_geo=fp_d
1100real(kind=fp_geo), PARAMETER :: fp_geo_miss=dmiss
1101
1104TYPE geo_coord
1105 PRIVATE
1106 INTEGER(kind=int_l) :: ilon, ilat
1107END TYPE geo_coord
1108
1109TYPE geo_coord_r
1110 PRIVATE
1111 REAL(kind=fp_geo) :: lon, lat
1112END TYPE geo_coord_r
1113
1114
1118TYPE geo_coordvect
1119 PRIVATE
1120 REAL(kind=fp_geo),POINTER :: ll(:,:) => null()
1121 INTEGER :: vsize, vtype
1122END TYPE geo_coordvect
1123
1125TYPE(geo_coord),PARAMETER :: geo_coord_miss=geo_coord(imiss,imiss)
1126
1127INTEGER, PARAMETER :: geo_coordvect_point = 1
1128INTEGER, PARAMETER :: geo_coordvect_arc = 3
1129INTEGER, PARAMETER :: geo_coordvect_polygon = 5
1130INTEGER, PARAMETER :: geo_coordvect_multipoint = 8
1131
1132
1136INTERFACE init
1137 MODULE PROCEDURE geo_coord_init, geo_coordvect_init
1138END INTERFACE
1139
1143INTERFACE delete
1144 MODULE PROCEDURE geo_coord_delete, geo_coordvect_delete
1145END INTERFACE
1146
1148INTERFACE getval
1149 MODULE PROCEDURE geo_coord_getval, geo_coordvect_getval
1150END INTERFACE
1151
1153INTERFACE OPERATOR (==)
1154 MODULE PROCEDURE geo_coord_eq
1155END INTERFACE
1156
1158INTERFACE OPERATOR (/=)
1159 MODULE PROCEDURE geo_coord_ne
1160END INTERFACE
1161
1162INTERFACE OPERATOR (>=)
1163 MODULE PROCEDURE geo_coord_ge
1164END INTERFACE
1165
1166INTERFACE OPERATOR (>)
1167 MODULE PROCEDURE geo_coord_gt
1168END INTERFACE
1169
1170INTERFACE OPERATOR (<=)
1171 MODULE PROCEDURE geo_coord_le
1172END INTERFACE
1173
1174INTERFACE OPERATOR (<)
1175 MODULE PROCEDURE geo_coord_lt
1176END INTERFACE
1177
1180INTERFACE import
1181 MODULE PROCEDURE geo_coordvect_import, geo_coordvect_importvect
1182END INTERFACE
1183
1186INTERFACE export
1187 MODULE PROCEDURE geo_coordvect_export, geo_coordvect_exportvect
1188END INTERFACE
1189
1192INTERFACE read_unit
1193 MODULE PROCEDURE geo_coord_read_unit, geo_coord_vect_read_unit
1194END INTERFACE
1195
1198INTERFACE write_unit
1199 MODULE PROCEDURE geo_coord_write_unit, geo_coord_vect_write_unit
1200END INTERFACE
1201
1203INTERFACE inside
1204 MODULE PROCEDURE geo_coord_inside, geo_coord_inside_rectang
1205END INTERFACE
1206
1208INTERFACE c_e
1209 MODULE PROCEDURE c_e_geo_coord
1210END INTERFACE
1211
1213INTERFACE to_char
1214 MODULE PROCEDURE to_char_geo_coord
1215END INTERFACE
1216
1218INTERFACE display
1219 MODULE PROCEDURE display_geo_coord
1220END INTERFACE
1221
1222CONTAINS
1223
1224
1225! ===================
1226! == geo_coord ==
1227! ===================
1234SUBROUTINE geo_coord_init(this, lon, lat, ilon, ilat)
1235TYPE(geo_coord) :: this
1236REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon
1237REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat
1238integer(kind=int_l), INTENT(IN), OPTIONAL :: ilon
1239integer(kind=int_l), INTENT(IN), OPTIONAL :: ilat
1240
1241real(kind=fp_geo) :: llon,llat
1242
1243CALL optio(ilon, this%ilon)
1244CALL optio(ilat, this%ilat)
1245
1246if (.not. c_e(this%ilon)) then
1247 CALL optio(lon, llon)
1248 if (c_e(llon)) this%ilon=nint(llon*1.d5)
1249end if
1250
1251if (.not. c_e(this%ilat)) then
1252 CALL optio(lat, llat)
1253 if (c_e(llat)) this%ilat=nint(llat*1.d5)
1254end if
1255
1256END SUBROUTINE geo_coord_init
1257
1259SUBROUTINE geo_coord_delete(this)
1260TYPE(geo_coord), INTENT(INOUT) :: this
1261
1262this%ilon = imiss
1263this%ilat = imiss
1264
1265END SUBROUTINE geo_coord_delete
1266
1271elemental SUBROUTINE geo_coord_getval(this, lon, lat,ilon,ilat)
1272TYPE(geo_coord),INTENT(IN) :: this
1273REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lon
1274REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lat
1275integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilon
1276integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilat
1277
1278IF (PRESENT(ilon)) ilon = getilon(this)
1279IF (PRESENT(ilat)) ilat = getilat(this)
1280
1281IF (PRESENT(lon)) lon = getlon(this)
1282IF (PRESENT(lat)) lat = getlat(this)
1283
1284END SUBROUTINE geo_coord_getval
1285
1286
1291elemental FUNCTION getilat(this)
1292TYPE(geo_coord),INTENT(IN) :: this
1293integer(kind=int_l) :: getilat
1294
1295getilat = this%ilat
1296
1297END FUNCTION getilat
1298
1303elemental FUNCTION getlat(this)
1304TYPE(geo_coord),INTENT(IN) :: this
1305real(kind=fp_geo) :: getlat
1306integer(kind=int_l) :: ilat
1307
1308ilat=getilat(this)
1309if (c_e(ilat)) then
1310 getlat = ilat*1.d-5
1311else
1312 getlat=fp_geo_miss
1313end if
1314
1315END FUNCTION getlat
1316
1321elemental FUNCTION getilon(this)
1322TYPE(geo_coord),INTENT(IN) :: this
1323integer(kind=int_l) :: getilon
1324
1325getilon = this%ilon
1326
1327END FUNCTION getilon
1328
1329
1334elemental FUNCTION getlon(this)
1335TYPE(geo_coord),INTENT(IN) :: this
1336real(kind=fp_geo) :: getlon
1337integer(kind=int_l) :: ilon
1338
1339ilon=getilon(this)
1340if (c_e(ilon)) then
1341 getlon = ilon*1.d-5
1342else
1343 getlon=fp_geo_miss
1344end if
1345
1346END FUNCTION getlon
1347
1348
1349elemental FUNCTION geo_coord_eq(this, that) RESULT(res)
1350TYPE(geo_coord),INTENT(IN) :: this, that
1351LOGICAL :: res
1352
1353res = (this%ilon == that%ilon .AND. this%ilat == that%ilat)
1354
1355END FUNCTION geo_coord_eq
1356
1357
1358elemental FUNCTION geo_coord_ne(this, that) RESULT(res)
1359TYPE(geo_coord),INTENT(IN) :: this, that
1360LOGICAL :: res
1361
1362res = .not. this == that
1363
1364END FUNCTION geo_coord_ne
1365
1368elemental FUNCTION geo_coord_gt(this, that) RESULT(res)
1369TYPE(geo_coord),INTENT(IN) :: this, that
1370LOGICAL :: res
1371
1372res = this%ilon > that%ilon
1373
1374if ( this%ilon == that%ilon ) then
1375 res = this%ilat > that%ilat
1376end if
1377
1378END FUNCTION geo_coord_gt
1379
1380
1383elemental FUNCTION geo_coord_ge(this, that) RESULT(res)
1384TYPE(geo_coord),INTENT(IN) :: this, that
1385LOGICAL :: res
1386
1387res = .not. this < that
1388
1389END FUNCTION geo_coord_ge
1390
1393elemental FUNCTION geo_coord_lt(this, that) RESULT(res)
1394TYPE(geo_coord),INTENT(IN) :: this, that
1395LOGICAL :: res
1396
1397res = this%ilon < that%ilon
1398
1399if ( this%ilon == that%ilon ) then
1400 res = this%ilat < that%ilat
1401end if
1402
1403
1404END FUNCTION geo_coord_lt
1405
1406
1409elemental FUNCTION geo_coord_le(this, that) RESULT(res)
1410TYPE(geo_coord),INTENT(IN) :: this, that
1411LOGICAL :: res
1412
1413res = .not. this > that
1414
1415END FUNCTION geo_coord_le
1416
1417
1420elemental FUNCTION geo_coord_ure(this, that) RESULT(res)
1421TYPE(geo_coord),INTENT(IN) :: this, that
1422LOGICAL :: res
1423
1424res = (this%ilon >= that%ilon .AND. this%ilat >= that%ilat)
1425
1426END FUNCTION geo_coord_ure
1427
1430elemental FUNCTION geo_coord_ur(this, that) RESULT(res)
1431TYPE(geo_coord),INTENT(IN) :: this, that
1432LOGICAL :: res
1433
1434res = (this%ilon > that%ilon .AND. this%ilat > that%ilat)
1435
1436END FUNCTION geo_coord_ur
1437
1438
1441elemental FUNCTION geo_coord_lle(this, that) RESULT(res)
1442TYPE(geo_coord),INTENT(IN) :: this, that
1443LOGICAL :: res
1444
1445res = (this%ilon <= that%ilon .AND. this%ilat <= that%ilat)
1446
1447END FUNCTION geo_coord_lle
1448
1451elemental FUNCTION geo_coord_ll(this, that) RESULT(res)
1452TYPE(geo_coord),INTENT(IN) :: this, that
1453LOGICAL :: res
1454
1455res = (this%ilon < that%ilon .AND. this%ilat < that%ilat)
1456
1457END FUNCTION geo_coord_ll
1458
1459
1465SUBROUTINE geo_coord_read_unit(this, unit)
1466TYPE(geo_coord),INTENT(out) :: this
1467INTEGER, INTENT(in) :: unit
1468
1469CALL geo_coord_vect_read_unit((/this/), unit)
1470
1471END SUBROUTINE geo_coord_read_unit
1472
1473
1479SUBROUTINE geo_coord_vect_read_unit(this, unit)
1480TYPE(geo_coord) :: this(:)
1481INTEGER, INTENT(in) :: unit
1482
1483CHARACTER(len=40) :: form
1484INTEGER :: i
1485
1486INQUIRE(unit, form=form)
1487IF (form == 'FORMATTED') THEN
1488 read(unit,*) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1489!TODO bug gfortran compiler !
1490!missing values are unredeable when formatted
1491ELSE
1492 READ(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1493ENDIF
1494
1495END SUBROUTINE geo_coord_vect_read_unit
1496
1497
1502SUBROUTINE geo_coord_write_unit(this, unit)
1503TYPE(geo_coord),INTENT(in) :: this
1504INTEGER, INTENT(in) :: unit
1505
1506CALL geo_coord_vect_write_unit((/this/), unit)
1507
1508END SUBROUTINE geo_coord_write_unit
1509
1510
1515SUBROUTINE geo_coord_vect_write_unit(this, unit)
1516TYPE(geo_coord),INTENT(in) :: this(:)
1517INTEGER, INTENT(in) :: unit
1518
1519CHARACTER(len=40) :: form
1520INTEGER :: i
1521
1522INQUIRE(unit, form=form)
1523IF (form == 'FORMATTED') THEN
1524 WRITE(unit,*) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1525!TODO bug gfortran compiler !
1526!missing values are unredeable when formatted
1527ELSE
1528 WRITE(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1529ENDIF
1530
1531END SUBROUTINE geo_coord_vect_write_unit
1532
1533
1536FUNCTION geo_coord_dist(this, that) RESULT(dist)
1538TYPE(geo_coord), INTENT (IN) :: this
1539TYPE(geo_coord), INTENT (IN) :: that
1540REAL(kind=fp_geo) :: dist
1541
1542REAL(kind=fp_geo) :: x,y
1543! Distanza approssimata, valida per piccoli angoli
1544
1545x = (getlon(this)-getlon(that))*cos(((getlat(this)+getlat(that))/2.)*degrad)
1546y = getlat(this)-getlat(that)
1547dist = sqrt(x**2 + y**2)*degrad*rearth
1548
1549END FUNCTION geo_coord_dist
1550
1551
1560FUNCTION geo_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1561TYPE(geo_coord),INTENT(IN) :: this
1562TYPE(geo_coord),INTENT(IN) :: coordmin
1563TYPE(geo_coord),INTENT(IN) :: coordmax
1564LOGICAL :: res
1565
1566res = (geo_coord_ure(this, coordmin) .AND. geo_coord_lle(this,coordmax))
1567
1568END FUNCTION geo_coord_inside_rectang
1569
1570
1571! ===================
1572! == geo_coordvect ==
1573! ===================
1582RECURSIVE SUBROUTINE geo_coordvect_init(this, lon, lat)
1583TYPE(geo_coordvect), INTENT(OUT) :: this
1584REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon(:)
1585REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat(:)
1586
1587IF (PRESENT(lon) .AND. PRESENT(lat)) THEN
1588 this%vsize = min(SIZE(lon), SIZE(lat))
1589 ALLOCATE(this%ll(this%vsize,2))
1590 this%ll(1:this%vsize,1) = lon(1:this%vsize)
1591 this%ll(1:this%vsize,2) = lat(1:this%vsize)
1592ELSE
1593 this%vsize = 0
1594 NULLIFY(this%ll)
1595ENDIF
1596this%vtype = 0 !?
1597
1598END SUBROUTINE geo_coordvect_init
1599
1600
1603SUBROUTINE geo_coordvect_delete(this)
1604TYPE(geo_coordvect), INTENT(INOUT) :: this
1605
1606IF (ASSOCIATED(this%ll)) DEALLOCATE(this%ll)
1607this%vsize = 0
1608this%vtype = 0
1609
1610END SUBROUTINE geo_coordvect_delete
1611
1612
1621SUBROUTINE geo_coordvect_getval(this, lon, lat)
1622TYPE(geo_coordvect),INTENT(IN) :: this
1623REAL(kind=fp_geo), OPTIONAL, POINTER :: lon(:)
1624REAL(kind=fp_geo), OPTIONAL, POINTER :: lat(:)
1625
1626IF (PRESENT(lon)) THEN
1627 IF (ASSOCIATED(this%ll)) THEN
1628 ALLOCATE(lon(this%vsize))
1629 lon(:) = this%ll(1:this%vsize,1)
1630 ENDIF
1631ENDIF
1632IF (PRESENT(lat)) THEN
1633 IF (ASSOCIATED(this%ll)) THEN
1634 ALLOCATE(lat(this%vsize))
1635 lat(:) = this%ll(1:this%vsize,2)
1636 ENDIF
1637ENDIF
1638
1639END SUBROUTINE geo_coordvect_getval
1640
1641
1642SUBROUTINE geo_coordvect_import(this, unitsim, &
1643#ifdef HAVE_SHAPELIB
1644 shphandle, &
1645#endif
1646 nshp)
1647TYPE(geo_coordvect), INTENT(OUT) :: this
1648INTEGER,OPTIONAL,INTENT(IN) :: unitsim
1649#ifdef HAVE_SHAPELIB
1650TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
1651#endif
1652INTEGER,OPTIONAL,INTENT(IN) :: nshp
1653
1654REAL(kind=fp_geo),ALLOCATABLE :: llon(:), llat(:)
1655REAL(kind=fp_geo) :: lv1,lv2,lv3,lv4
1656INTEGER :: i, lvsize
1657CHARACTER(len=40) :: lname
1658#ifdef HAVE_SHAPELIB
1659TYPE(shpobject) :: shpobj
1660#endif
1661
1662IF (PRESENT(unitsim)) THEN
1663 ! Leggo l'intestazione
1664 READ(unitsim,*,END=10)lvsize,lv1,lv2,lv3,lv4,lname
1665 ALLOCATE(llon(lvsize+1), llat(lvsize+1))
1666 ! Leggo il poligono
1667 READ(unitsim,*)(llon(i),llat(i), i=1,lvsize)
1668 ! Lo chiudo se necessario
1669 IF (llon(1) /= llon(lvsize) .OR. llat(1) /= llat(lvsize)) THEN
1670 lvsize = lvsize + 1
1671 llon(lvsize) = llon(1)
1672 llat(lvsize) = llat(1)
1673 ENDIF
1674 ! Lo inserisco nel mio oggetto
1675 CALL init(this, lon=llon(1:lvsize), lat=llat(1:lvsize))
1676 this%vtype = geo_coordvect_polygon ! Sempre un poligono
1677
1678 DEALLOCATE(llon, llat)
1679 RETURN
168010 CALL raise_error('nella lettura del file '//trim(to_char(unitsim)))
1681 DEALLOCATE(llon, llat) ! End of file, ritorno un oggetto non assegnato
1682#ifdef HAVE_SHAPELIB
1683ELSE IF (PRESENT(shphandle) .AND. PRESENT(nshp)) THEN
1684 ! Leggo l'oggetto shape
1685 shpobj = shpreadobject(shphandle, nshp)
1686 IF (.NOT.shpisnull(shpobj)) THEN
1687 ! Lo inserisco nel mio oggetto
1688 CALL init(this, lon=real(shpobj%padfx,kind=fp_geo), &
1689 lat=real(shpobj%padfy,kind=fp_geo))
1690 this%vtype = shpobj%nshptype
1691 CALL shpdestroyobject(shpobj)
1692 ELSE
1693 CALL init(this)
1694 ENDIF
1695#endif
1696ENDIF
1697
1698END SUBROUTINE geo_coordvect_import
1699
1700
1701SUBROUTINE geo_coordvect_export(this, unitsim, &
1702#ifdef HAVE_SHAPELIB
1703 shphandle, &
1704#endif
1705 nshp)
1706TYPE(geo_coordvect), INTENT(INOUT) :: this
1707INTEGER,OPTIONAL,INTENT(IN) :: unitsim
1708#ifdef HAVE_SHAPELIB
1709TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
1710#endif
1711INTEGER,OPTIONAL,INTENT(IN) :: nshp
1712
1713INTEGER :: i, lnshp
1714#ifdef HAVE_SHAPELIB
1715TYPE(shpobject) :: shpobj
1716#endif
1717
1718IF (PRESENT(unitsim)) THEN
1719 IF (this%vsize > 0) THEN
1720 ! Scrivo l'intestazione
1721 WRITE(unitsim,*)SIZE(this%ll,1),-1.,5000.,-0.1,1.1,'Area'
1722 ! Scrivo il poligono
1723 WRITE(unitsim,*)(this%ll(i,1:2), i=1,this%vsize)
1724 ELSE
1725 CALL raise_warning('oggetto geo_coordvect vuoto, non scrivo niente in '// &
1726 trim(to_char(unitsim)))
1727 ENDIF
1728#ifdef HAVE_SHAPELIB
1729ELSE IF (PRESENT(shphandle)) THEN
1730 IF (PRESENT(nshp)) THEN
1731 lnshp = nshp
1732 ELSE
1733 lnshp = -1 ! -1 = append
1734 ENDIF
1735 ! Creo l'oggetto shape inizializzandolo con il mio oggetto
1736 shpobj = shpcreatesimpleobject(this%vtype, this%vsize, &
1737 REAL(this%ll(1:this%vsize,1),kind=fp_d), &
1738 REAL(this%ll(1:this%vsize,2),kind=fp_d))
1739 IF (.NOT.shpisnull(shpobj)) THEN
1740 ! Lo scrivo nello shapefile
1741 i=shpwriteobject(shphandle, lnshp, shpobj)
1742 CALL shpdestroyobject(shpobj)
1743 ENDIF
1744#endif
1745ENDIF
1746
1747END SUBROUTINE geo_coordvect_export
1748
1767SUBROUTINE geo_coordvect_importvect(this, shpfilesim, shpfile)
1768TYPE(geo_coordvect),POINTER :: this(:)
1769CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim
1770CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile
1771
1772REAL(kind=fp_geo) :: inu
1773REAL(kind=fp_d) :: minb(4), maxb(4)
1774INTEGER :: i, u, ns, lvsize, shptype, dbfnf, dbfnr
1775CHARACTER(len=40) :: lname
1776#ifdef HAVE_SHAPELIB
1777TYPE(shpfileobject) :: shphandle
1778#endif
1779
1780NULLIFY(this)
1781
1782IF (PRESENT(shpfilesim)) THEN
1783 u = getunit()
1784 OPEN(u, file=shpfilesim, status='old', err=30)
1785 ns = 0 ! Conto il numero di shape contenute
1786 DO WHILE(.true.)
1787 READ(u,*,END=10,ERR=20)lvsize,inu,inu,inu,inu,lname
1788 READ(u,*,END=20,ERR=20)(inu,inu, i=1,lvsize)
1789 ns = ns + 1
1790 ENDDO
179110 CONTINUE
1792 IF (ns > 0) THEN ! Alloco e leggo il mio oggetto
1793 ALLOCATE(this(ns))
1794 rewind(u)
1795 DO i = 1, ns
1796 CALL import(this(i), unitsim=u)
1797 ENDDO
1798 ENDIF
179920 CONTINUE
1800 CLOSE(u)
1801 IF (.NOT.ASSOCIATED(this)) THEN
1802 CALL raise_warning('file '//trim(shpfilesim)//' vuoto o corrotto')
1803 ENDIF
1804 RETURN
180530 CONTINUE
1806 CALL raise_error('Impossibile aprire il file '//trim(shpfile))
1807 RETURN
1808
1809ELSE IF (PRESENT(shpfile)) THEN
1810#ifdef HAVE_SHAPELIB
1811 shphandle = shpopen(trim(shpfile), 'rb')
1812 IF (shpfileisnull(shphandle)) THEN
1813 CALL raise_error('Impossibile aprire lo shapefile '//trim(shpfile))
1814 RETURN
1815 ENDIF
1816 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr) ! Ottengo le info sul file
1817 IF (ns > 0) THEN ! Alloco e leggo il mio oggetto
1818 ALLOCATE(this(ns))
1819 this(:)%vtype = shptype
1820 DO i = 1, ns
1821 CALL import(this(i), shphandle=shphandle, nshp=i-1)
1822 ENDDO
1823 ENDIF
1824 CALL shpclose(shphandle)
1825 RETURN
1826#endif
1827ENDIF
1828
1829END SUBROUTINE geo_coordvect_importvect
1830
1831
1834SUBROUTINE geo_coordvect_exportvect(this, shpfilesim, shpfile, append)
1835TYPE(geo_coordvect) :: this(:)
1836CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim
1837CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile
1843LOGICAL,INTENT(in),OPTIONAL :: append
1844
1845REAL(kind=fp_d) :: minb(4), maxb(4)
1846INTEGER :: i, u, ns, shptype, dbfnf, dbfnr
1847LOGICAL :: lappend
1848#ifdef HAVE_SHAPELIB
1849TYPE(shpfileobject) :: shphandle
1850#endif
1851
1852IF (PRESENT(append)) THEN
1853 lappend = append
1854ELSE
1855 lappend = .false.
1856ENDIF
1857IF (PRESENT(shpfilesim)) THEN
1858 u = getunit()
1859 IF (lappend) THEN
1860 OPEN(u, file=shpfilesim, status='unknown', position='append', err=30)
1861 ELSE
1862 OPEN(u, file=shpfilesim, status='unknown', err=30)
1863 ENDIF
1864 DO i = 1, SIZE(this)
1865 CALL export(this(i), unitsim=u)
1866 ENDDO
1867 CLOSE(u)
1868 RETURN
186930 CONTINUE
1870 CALL raise_error('Impossibile aprire il file '//trim(shpfile))
1871 RETURN
1872ELSE IF (PRESENT(shpfile)) THEN
1873#ifdef HAVE_SHAPELIB
1874 IF (lappend) THEN
1875 shphandle = shpopen(trim(shpfile), 'r+b')
1876 IF (shpfileisnull(shphandle)) THEN ! shpopen funziona solo su file esistenti
1877 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
1878 ENDIF
1879 ELSE
1880 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
1881 ENDIF
1882 IF (shpfileisnull(shphandle)) THEN
1883 CALL raise_error('Impossibile aprire lo shapefile '//trim(shpfile))
1884 RETURN
1885 ENDIF
1886 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr) ! Ottengo le info sul file
1887 DO i = 1, SIZE(this)
1888 IF (i > ns .OR. lappend) THEN ! Append shape
1889 CALL export(this(i), shphandle=shphandle)
1890 ELSE ! Overwrite shape
1891 CALL export(this(i), shphandle=shphandle, nshp=i-1)
1892 ENDIF
1893 ENDDO
1894 CALL shpclose(shphandle)
1895 RETURN
1896#endif
1897ENDIF
1898
1899END SUBROUTINE geo_coordvect_exportvect
1900
1901
1910FUNCTION geo_coord_inside(this, poly) RESULT(inside)
1911TYPE(geo_coord), INTENT(IN) :: this
1912TYPE(geo_coordvect), INTENT(IN) :: poly
1913LOGICAL :: inside
1914
1915INTEGER :: i, j, starti
1916
1917inside = .false.
1918IF (all(poly%ll(1,:) == poly%ll(poly%vsize,:))) THEN ! Poligono chiuso
1919 starti = 2
1920 j = 1
1921ELSE ! Poligono non chiuso
1922 starti = 1
1923 j = poly%vsize
1924ENDIF
1925DO i = starti, poly%vsize
1926 IF ((poly%ll(i,2) <= getlat(this) .AND. &
1927 getlat(this) < poly%ll(j,2)) .OR. &
1928 (poly%ll(j,2) <= getlat(this) .AND. &
1929 getlat(this) < poly%ll(i,2))) THEN
1930 IF (getlon(this) < (poly%ll(j,1) - poly%ll(i,1)) * &
1931 (getlat(this) - poly%ll(i,2)) / &
1932 (poly%ll(j,2) - poly%ll(i,2)) + poly%ll(i,1)) THEN
1933 inside = .NOT. inside
1934 ENDIF
1935 ENDIF
1936 j = i
1937ENDDO
1938
1939END FUNCTION geo_coord_inside
1940
1941
1942ELEMENTAL FUNCTION c_e_geo_coord(this) result (res)
1943TYPE(geo_coord),INTENT(in) :: this
1944LOGICAL :: res
1945
1946res = .not. this == geo_coord_miss
1947
1948end FUNCTION c_e_geo_coord
1949
1950
1951character(len=80) function to_char_geo_coord(this)
1952TYPE(geo_coord),INTENT(in) :: this
1953
1954to_char_geo_coord = "GEO_COORD: Lon="// &
1955 trim(to_char(getlon(this),miss="Missing lon",form="(f11.5)"))//&
1956 " Lat="// &
1957 trim(to_char(getlat(this),miss="Missing lat",form="(f11.5)"))
1958
1959end function to_char_geo_coord
1960
1961
1962subroutine display_geo_coord(this)
1963TYPE(geo_coord),INTENT(in) :: this
1964
1965print*,trim(to_char(this))
1966
1967end subroutine display_geo_coord
1968
1969
1970END 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:245
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.