libsim Versione 7.1.11
|
◆ geo_coordvect_exportvect()
Esporta un vettore di oggetti geo_coordvect su un file in formato testo o in formato shapefile.
Definizione alla linea 979 del file geo_coord_class.F90. 980! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
981! authors:
982! Davide Cesari <dcesari@arpa.emr.it>
983! Paolo Patruno <ppatruno@arpa.emr.it>
984
985! This program is free software; you can redistribute it and/or
986! modify it under the terms of the GNU General Public License as
987! published by the Free Software Foundation; either version 2 of
988! the License, or (at your option) any later version.
989
990! This program is distributed in the hope that it will be useful,
991! but WITHOUT ANY WARRANTY; without even the implied warranty of
992! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
993! GNU General Public License for more details.
994
995! You should have received a copy of the GNU General Public License
996! along with this program. If not, see <http://www.gnu.org/licenses/>.
997#include "config.h"
998
1012!USE doubleprecision_phys_const
1014#ifdef HAVE_SHAPELIB
1015USE shapelib
1016!, ONLY: shpobject, shpreadobject, shpdestroyobject, shpwriteobject, &
1017! shpcreatesimpleobject
1018#endif
1019IMPLICIT NONE
1020
1021
1029INTEGER, PARAMETER :: fp_geo=fp_d
1030real(kind=fp_geo), PARAMETER :: fp_geo_miss=dmiss
1031
1035 PRIVATE
1036 INTEGER(kind=int_l) :: ilon, ilat
1038
1039TYPE geo_coord_r
1040 PRIVATE
1041 REAL(kind=fp_geo) :: lon, lat
1042END TYPE geo_coord_r
1043
1044
1049 PRIVATE
1050 REAL(kind=fp_geo),POINTER :: ll(:,:) => null()
1051 INTEGER :: vsize, vtype
1053
1056
1057INTEGER, PARAMETER :: geo_coordvect_point = 1
1058INTEGER, PARAMETER :: geo_coordvect_arc = 3
1059INTEGER, PARAMETER :: geo_coordvect_polygon = 5
1060INTEGER, PARAMETER :: geo_coordvect_multipoint = 8
1061
1062
1067 MODULE PROCEDURE geo_coord_init, geo_coordvect_init
1068END INTERFACE
1069
1074 MODULE PROCEDURE geo_coord_delete, geo_coordvect_delete
1075END INTERFACE
1076
1079 MODULE PROCEDURE geo_coord_getval, geo_coordvect_getval
1080END INTERFACE
1081
1083INTERFACE OPERATOR (==)
1084 MODULE PROCEDURE geo_coord_eq
1085END INTERFACE
1086
1088INTERFACE OPERATOR (/=)
1089 MODULE PROCEDURE geo_coord_ne
1090END INTERFACE
1091
1092INTERFACE OPERATOR (>=)
1093 MODULE PROCEDURE geo_coord_ge
1094END INTERFACE
1095
1096INTERFACE OPERATOR (>)
1097 MODULE PROCEDURE geo_coord_gt
1098END INTERFACE
1099
1100INTERFACE OPERATOR (<=)
1101 MODULE PROCEDURE geo_coord_le
1102END INTERFACE
1103
1104INTERFACE OPERATOR (<)
1105 MODULE PROCEDURE geo_coord_lt
1106END INTERFACE
1107
1110INTERFACE import
1111 MODULE PROCEDURE geo_coordvect_import, geo_coordvect_importvect
1112END INTERFACE
1113
1117 MODULE PROCEDURE geo_coordvect_export, geo_coordvect_exportvect
1118END INTERFACE
1119
1123 MODULE PROCEDURE geo_coord_read_unit, geo_coord_vect_read_unit
1124END INTERFACE
1125
1129 MODULE PROCEDURE geo_coord_write_unit, geo_coord_vect_write_unit
1130END INTERFACE
1131
1134 MODULE PROCEDURE geo_coord_inside, geo_coord_inside_rectang
1135END INTERFACE
1136
1139 MODULE PROCEDURE c_e_geo_coord
1140END INTERFACE
1141
1144 MODULE PROCEDURE to_char_geo_coord
1145END INTERFACE
1146
1149 MODULE PROCEDURE display_geo_coord
1150END INTERFACE
1151
1152CONTAINS
1153
1154
1155! ===================
1156! == geo_coord ==
1157! ===================
1164SUBROUTINE geo_coord_init(this, lon, lat, ilon, ilat)
1165TYPE(geo_coord) :: this
1166REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon
1167REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat
1168integer(kind=int_l), INTENT(IN), OPTIONAL :: ilon
1169integer(kind=int_l), INTENT(IN), OPTIONAL :: ilat
1170
1171real(kind=fp_geo) :: llon,llat
1172
1173CALL optio(ilon, this%ilon)
1174CALL optio(ilat, this%ilat)
1175
1177 CALL optio(lon, llon)
1179end if
1180
1182 CALL optio(lat, llat)
1184end if
1185
1186END SUBROUTINE geo_coord_init
1187
1189SUBROUTINE geo_coord_delete(this)
1190TYPE(geo_coord), INTENT(INOUT) :: this
1191
1192this%ilon = imiss
1193this%ilat = imiss
1194
1195END SUBROUTINE geo_coord_delete
1196
1201elemental SUBROUTINE geo_coord_getval(this, lon, lat,ilon,ilat)
1202TYPE(geo_coord),INTENT(IN) :: this
1203REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lon
1204REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lat
1205integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilon
1206integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilat
1207
1208IF (PRESENT(ilon)) ilon = getilon(this)
1209IF (PRESENT(ilat)) ilat = getilat(this)
1210
1211IF (PRESENT(lon)) lon = getlon(this)
1212IF (PRESENT(lat)) lat = getlat(this)
1213
1214END SUBROUTINE geo_coord_getval
1215
1216
1221elemental FUNCTION getilat(this)
1222TYPE(geo_coord),INTENT(IN) :: this
1223integer(kind=int_l) :: getilat
1224
1225getilat = this%ilat
1226
1227END FUNCTION getilat
1228
1233elemental FUNCTION getlat(this)
1234TYPE(geo_coord),INTENT(IN) :: this
1235real(kind=fp_geo) :: getlat
1236integer(kind=int_l) :: ilat
1237
1238ilat=getilat(this)
1240 getlat = ilat*1.d-5
1241else
1242 getlat=fp_geo_miss
1243end if
1244
1245END FUNCTION getlat
1246
1251elemental FUNCTION getilon(this)
1252TYPE(geo_coord),INTENT(IN) :: this
1253integer(kind=int_l) :: getilon
1254
1255getilon = this%ilon
1256
1257END FUNCTION getilon
1258
1259
1264elemental FUNCTION getlon(this)
1265TYPE(geo_coord),INTENT(IN) :: this
1266real(kind=fp_geo) :: getlon
1267integer(kind=int_l) :: ilon
1268
1269ilon=getilon(this)
1271 getlon = ilon*1.d-5
1272else
1273 getlon=fp_geo_miss
1274end if
1275
1276END FUNCTION getlon
1277
1278
1279elemental FUNCTION geo_coord_eq(this, that) RESULT(res)
1280TYPE(geo_coord),INTENT(IN) :: this, that
1281LOGICAL :: res
1282
1283res = (this%ilon == that%ilon .AND. this%ilat == that%ilat)
1284
1285END FUNCTION geo_coord_eq
1286
1287
1288elemental FUNCTION geo_coord_ne(this, that) RESULT(res)
1289TYPE(geo_coord),INTENT(IN) :: this, that
1290LOGICAL :: res
1291
1292res = .not. this == that
1293
1294END FUNCTION geo_coord_ne
1295
1298elemental FUNCTION geo_coord_gt(this, that) RESULT(res)
1299TYPE(geo_coord),INTENT(IN) :: this, that
1300LOGICAL :: res
1301
1302res = this%ilon > that%ilon
1303
1304if ( this%ilon == that%ilon ) then
1305 res = this%ilat > that%ilat
1306end if
1307
1308END FUNCTION geo_coord_gt
1309
1310
1313elemental FUNCTION geo_coord_ge(this, that) RESULT(res)
1314TYPE(geo_coord),INTENT(IN) :: this, that
1315LOGICAL :: res
1316
1317res = .not. this < that
1318
1319END FUNCTION geo_coord_ge
1320
1323elemental FUNCTION geo_coord_lt(this, that) RESULT(res)
1324TYPE(geo_coord),INTENT(IN) :: this, that
1325LOGICAL :: res
1326
1327res = this%ilon < that%ilon
1328
1329if ( this%ilon == that%ilon ) then
1330 res = this%ilat < that%ilat
1331end if
1332
1333
1334END FUNCTION geo_coord_lt
1335
1336
1339elemental FUNCTION geo_coord_le(this, that) RESULT(res)
1340TYPE(geo_coord),INTENT(IN) :: this, that
1341LOGICAL :: res
1342
1343res = .not. this > that
1344
1345END FUNCTION geo_coord_le
1346
1347
1350elemental FUNCTION geo_coord_ure(this, that) RESULT(res)
1351TYPE(geo_coord),INTENT(IN) :: this, that
1352LOGICAL :: res
1353
1354res = (this%ilon >= that%ilon .AND. this%ilat >= that%ilat)
1355
1356END FUNCTION geo_coord_ure
1357
1360elemental FUNCTION geo_coord_ur(this, that) RESULT(res)
1361TYPE(geo_coord),INTENT(IN) :: this, that
1362LOGICAL :: res
1363
1364res = (this%ilon > that%ilon .AND. this%ilat > that%ilat)
1365
1366END FUNCTION geo_coord_ur
1367
1368
1371elemental FUNCTION geo_coord_lle(this, that) RESULT(res)
1372TYPE(geo_coord),INTENT(IN) :: this, that
1373LOGICAL :: res
1374
1375res = (this%ilon <= that%ilon .AND. this%ilat <= that%ilat)
1376
1377END FUNCTION geo_coord_lle
1378
1381elemental FUNCTION geo_coord_ll(this, that) RESULT(res)
1382TYPE(geo_coord),INTENT(IN) :: this, that
1383LOGICAL :: res
1384
1385res = (this%ilon < that%ilon .AND. this%ilat < that%ilat)
1386
1387END FUNCTION geo_coord_ll
1388
1389
1395SUBROUTINE geo_coord_read_unit(this, unit)
1396TYPE(geo_coord),INTENT(out) :: this
1397INTEGER, INTENT(in) :: unit
1398
1399CALL geo_coord_vect_read_unit((/this/), unit)
1400
1401END SUBROUTINE geo_coord_read_unit
1402
1403
1409SUBROUTINE geo_coord_vect_read_unit(this, unit)
1410TYPE(geo_coord) :: this(:)
1411INTEGER, INTENT(in) :: unit
1412
1413CHARACTER(len=40) :: form
1414INTEGER :: i
1415
1416INQUIRE(unit, form=form)
1417IF (form == 'FORMATTED') THEN
1418 read(unit,*) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1419!TODO bug gfortran compiler !
1420!missing values are unredeable when formatted
1421ELSE
1422 READ(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1423ENDIF
1424
1425END SUBROUTINE geo_coord_vect_read_unit
1426
1427
1432SUBROUTINE geo_coord_write_unit(this, unit)
1433TYPE(geo_coord),INTENT(in) :: this
1434INTEGER, INTENT(in) :: unit
1435
1436CALL geo_coord_vect_write_unit((/this/), unit)
1437
1438END SUBROUTINE geo_coord_write_unit
1439
1440
1445SUBROUTINE geo_coord_vect_write_unit(this, unit)
1446TYPE(geo_coord),INTENT(in) :: this(:)
1447INTEGER, INTENT(in) :: unit
1448
1449CHARACTER(len=40) :: form
1450INTEGER :: i
1451
1452INQUIRE(unit, form=form)
1453IF (form == 'FORMATTED') THEN
1454 WRITE(unit,*) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1455!TODO bug gfortran compiler !
1456!missing values are unredeable when formatted
1457ELSE
1458 WRITE(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1459ENDIF
1460
1461END SUBROUTINE geo_coord_vect_write_unit
1462
1463
1466FUNCTION geo_coord_dist(this, that) RESULT(dist)
1468TYPE(geo_coord), INTENT (IN) :: this
1469TYPE(geo_coord), INTENT (IN) :: that
1470REAL(kind=fp_geo) :: dist
1471
1472REAL(kind=fp_geo) :: x,y
1473! Distanza approssimata, valida per piccoli angoli
1474
1475x = (getlon(this)-getlon(that))*cos(((getlat(this)+getlat(that))/2.)*degrad)
1476y = getlat(this)-getlat(that)
1477dist = sqrt(x**2 + y**2)*degrad*rearth
1478
1479END FUNCTION geo_coord_dist
1480
1481
1490FUNCTION geo_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1491TYPE(geo_coord),INTENT(IN) :: this
1492TYPE(geo_coord),INTENT(IN) :: coordmin
1493TYPE(geo_coord),INTENT(IN) :: coordmax
1494LOGICAL :: res
1495
1496res = (geo_coord_ure(this, coordmin) .AND. geo_coord_lle(this,coordmax))
1497
1498END FUNCTION geo_coord_inside_rectang
1499
1500
1501! ===================
1502! == geo_coordvect ==
1503! ===================
1512RECURSIVE SUBROUTINE geo_coordvect_init(this, lon, lat)
1513TYPE(geo_coordvect), INTENT(OUT) :: this
1514REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon(:)
1515REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat(:)
1516
1517IF (PRESENT(lon) .AND. PRESENT(lat)) THEN
1518 this%vsize = min(SIZE(lon), SIZE(lat))
1519 ALLOCATE(this%ll(this%vsize,2))
1520 this%ll(1:this%vsize,1) = lon(1:this%vsize)
1521 this%ll(1:this%vsize,2) = lat(1:this%vsize)
1522ELSE
1523 this%vsize = 0
1524 NULLIFY(this%ll)
1525ENDIF
1526this%vtype = 0 !?
1527
1528END SUBROUTINE geo_coordvect_init
1529
1530
1533SUBROUTINE geo_coordvect_delete(this)
1534TYPE(geo_coordvect), INTENT(INOUT) :: this
1535
1536IF (ASSOCIATED(this%ll)) DEALLOCATE(this%ll)
1537this%vsize = 0
1538this%vtype = 0
1539
1540END SUBROUTINE geo_coordvect_delete
1541
1542
1551SUBROUTINE geo_coordvect_getval(this, lon, lat)
1552TYPE(geo_coordvect),INTENT(IN) :: this
1553REAL(kind=fp_geo), OPTIONAL, POINTER :: lon(:)
1554REAL(kind=fp_geo), OPTIONAL, POINTER :: lat(:)
1555
1556IF (PRESENT(lon)) THEN
1557 IF (ASSOCIATED(this%ll)) THEN
1558 ALLOCATE(lon(this%vsize))
1559 lon(:) = this%ll(1:this%vsize,1)
1560 ENDIF
1561ENDIF
1562IF (PRESENT(lat)) THEN
1563 IF (ASSOCIATED(this%ll)) THEN
1564 ALLOCATE(lat(this%vsize))
1565 lat(:) = this%ll(1:this%vsize,2)
1566 ENDIF
1567ENDIF
1568
1569END SUBROUTINE geo_coordvect_getval
1570
1571
1572SUBROUTINE geo_coordvect_import(this, unitsim, &
1573#ifdef HAVE_SHAPELIB
1574 shphandle, &
1575#endif
1576 nshp)
1577TYPE(geo_coordvect), INTENT(OUT) :: this
1578INTEGER,OPTIONAL,INTENT(IN) :: unitsim
1579#ifdef HAVE_SHAPELIB
1580TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
1581#endif
1582INTEGER,OPTIONAL,INTENT(IN) :: nshp
1583
1584REAL(kind=fp_geo),ALLOCATABLE :: llon(:), llat(:)
1585REAL(kind=fp_geo) :: lv1,lv2,lv3,lv4
1586INTEGER :: i, lvsize
1587CHARACTER(len=40) :: lname
1588#ifdef HAVE_SHAPELIB
1589TYPE(shpobject) :: shpobj
1590#endif
1591
1592IF (PRESENT(unitsim)) THEN
1593 ! Leggo l'intestazione
1594 READ(unitsim,*,END=10)lvsize,lv1,lv2,lv3,lv4,lname
1595 ALLOCATE(llon(lvsize+1), llat(lvsize+1))
1596 ! Leggo il poligono
1597 READ(unitsim,*)(llon(i),llat(i), i=1,lvsize)
1598 ! Lo chiudo se necessario
1599 IF (llon(1) /= llon(lvsize) .OR. llat(1) /= llat(lvsize)) THEN
1600 lvsize = lvsize + 1
1601 llon(lvsize) = llon(1)
1602 llat(lvsize) = llat(1)
1603 ENDIF
1604 ! Lo inserisco nel mio oggetto
1606 this%vtype = geo_coordvect_polygon ! Sempre un poligono
1607
1608 DEALLOCATE(llon, llat)
1609 RETURN
1611 DEALLOCATE(llon, llat) ! End of file, ritorno un oggetto non assegnato
1612#ifdef HAVE_SHAPELIB
1613ELSE IF (PRESENT(shphandle) .AND. PRESENT(nshp)) THEN
1614 ! Leggo l'oggetto shape
1615 shpobj = shpreadobject(shphandle, nshp)
1616 IF (.NOT.shpisnull(shpobj)) THEN
1617 ! Lo inserisco nel mio oggetto
1619 lat=real(shpobj%padfy,kind=fp_geo))
1620 this%vtype = shpobj%nshptype
1621 CALL shpdestroyobject(shpobj)
1622 ELSE
1624 ENDIF
1625#endif
1626ENDIF
1627
1628END SUBROUTINE geo_coordvect_import
1629
1630
1631SUBROUTINE geo_coordvect_export(this, unitsim, &
1632#ifdef HAVE_SHAPELIB
1633 shphandle, &
1634#endif
1635 nshp)
1636TYPE(geo_coordvect), INTENT(INOUT) :: this
1637INTEGER,OPTIONAL,INTENT(IN) :: unitsim
1638#ifdef HAVE_SHAPELIB
1639TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
1640#endif
1641INTEGER,OPTIONAL,INTENT(IN) :: nshp
1642
1643INTEGER :: i, lnshp
1644#ifdef HAVE_SHAPELIB
1645TYPE(shpobject) :: shpobj
1646#endif
1647
1648IF (PRESENT(unitsim)) THEN
1649 IF (this%vsize > 0) THEN
1650 ! Scrivo l'intestazione
1651 WRITE(unitsim,*)SIZE(this%ll,1),-1.,5000.,-0.1,1.1,'Area'
1652 ! Scrivo il poligono
1653 WRITE(unitsim,*)(this%ll(i,1:2), i=1,this%vsize)
1654 ELSE
1655 CALL raise_warning('oggetto geo_coordvect vuoto, non scrivo niente in '// &
1656 trim(to_char(unitsim)))
1657 ENDIF
1658#ifdef HAVE_SHAPELIB
1659ELSE IF (PRESENT(shphandle)) THEN
1660 IF (PRESENT(nshp)) THEN
1661 lnshp = nshp
1662 ELSE
1663 lnshp = -1 ! -1 = append
1664 ENDIF
1665 ! Creo l'oggetto shape inizializzandolo con il mio oggetto
1666 shpobj = shpcreatesimpleobject(this%vtype, this%vsize, &
1667 REAL(this%ll(1:this%vsize,1),kind=fp_d), &
1668 REAL(this%ll(1:this%vsize,2),kind=fp_d))
1669 IF (.NOT.shpisnull(shpobj)) THEN
1670 ! Lo scrivo nello shapefile
1671 i=shpwriteobject(shphandle, lnshp, shpobj)
1672 CALL shpdestroyobject(shpobj)
1673 ENDIF
1674#endif
1675ENDIF
1676
1677END SUBROUTINE geo_coordvect_export
1678
1697SUBROUTINE geo_coordvect_importvect(this, shpfilesim, shpfile)
1698TYPE(geo_coordvect),POINTER :: this(:)
1699CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim
1700CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile
1701
1702REAL(kind=fp_geo) :: inu
1703REAL(kind=fp_d) :: minb(4), maxb(4)
1704INTEGER :: i, u, ns, lvsize, shptype, dbfnf, dbfnr
1705CHARACTER(len=40) :: lname
1706#ifdef HAVE_SHAPELIB
1707TYPE(shpfileobject) :: shphandle
1708#endif
1709
1710NULLIFY(this)
1711
1712IF (PRESENT(shpfilesim)) THEN
1713 u = getunit()
1714 OPEN(u, file=shpfilesim, status='old', err=30)
1715 ns = 0 ! Conto il numero di shape contenute
1716 DO WHILE(.true.)
1717 READ(u,*,END=10,ERR=20)lvsize,inu,inu,inu,inu,lname
1718 READ(u,*,END=20,ERR=20)(inu,inu, i=1,lvsize)
1719 ns = ns + 1
1720 ENDDO
172110 CONTINUE
1722 IF (ns > 0) THEN ! Alloco e leggo il mio oggetto
1723 ALLOCATE(this(ns))
1724 rewind(u)
1725 DO i = 1, ns
1727 ENDDO
1728 ENDIF
172920 CONTINUE
1730 CLOSE(u)
1731 IF (.NOT.ASSOCIATED(this)) THEN
1732 CALL raise_warning('file '//trim(shpfilesim)//' vuoto o corrotto')
1733 ENDIF
1734 RETURN
173530 CONTINUE
1736 CALL raise_error('Impossibile aprire il file '//trim(shpfile))
1737 RETURN
1738
1739ELSE IF (PRESENT(shpfile)) THEN
1740#ifdef HAVE_SHAPELIB
1741 shphandle = shpopen(trim(shpfile), 'rb')
1742 IF (shpfileisnull(shphandle)) THEN
1743 CALL raise_error('Impossibile aprire lo shapefile '//trim(shpfile))
1744 RETURN
1745 ENDIF
1746 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr) ! Ottengo le info sul file
1747 IF (ns > 0) THEN ! Alloco e leggo il mio oggetto
1748 ALLOCATE(this(ns))
1749 this(:)%vtype = shptype
1750 DO i = 1, ns
1752 ENDDO
1753 ENDIF
1754 CALL shpclose(shphandle)
1755 RETURN
1756#endif
1757ENDIF
1758
1759END SUBROUTINE geo_coordvect_importvect
1760
1761
1764SUBROUTINE geo_coordvect_exportvect(this, shpfilesim, shpfile, append)
1765TYPE(geo_coordvect) :: this(:)
1766CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim
1767CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile
1773LOGICAL,INTENT(in),OPTIONAL :: append
1774
1775REAL(kind=fp_d) :: minb(4), maxb(4)
1776INTEGER :: i, u, ns, shptype, dbfnf, dbfnr
1777LOGICAL :: lappend
1778#ifdef HAVE_SHAPELIB
1779TYPE(shpfileobject) :: shphandle
1780#endif
1781
1782IF (PRESENT(append)) THEN
1783 lappend = append
1784ELSE
1785 lappend = .false.
1786ENDIF
1787IF (PRESENT(shpfilesim)) THEN
1788 u = getunit()
1789 IF (lappend) THEN
1790 OPEN(u, file=shpfilesim, status='unknown', position='append', err=30)
1791 ELSE
1792 OPEN(u, file=shpfilesim, status='unknown', err=30)
1793 ENDIF
1794 DO i = 1, SIZE(this)
1796 ENDDO
1797 CLOSE(u)
1798 RETURN
179930 CONTINUE
1800 CALL raise_error('Impossibile aprire il file '//trim(shpfile))
1801 RETURN
1802ELSE IF (PRESENT(shpfile)) THEN
1803#ifdef HAVE_SHAPELIB
1804 IF (lappend) THEN
1805 shphandle = shpopen(trim(shpfile), 'r+b')
1806 IF (shpfileisnull(shphandle)) THEN ! shpopen funziona solo su file esistenti
1807 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
1808 ENDIF
1809 ELSE
1810 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
1811 ENDIF
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 DO i = 1, SIZE(this)
1818 IF (i > ns .OR. lappend) THEN ! Append shape
1820 ELSE ! Overwrite shape
1822 ENDIF
1823 ENDDO
1824 CALL shpclose(shphandle)
1825 RETURN
1826#endif
1827ENDIF
1828
1829END SUBROUTINE geo_coordvect_exportvect
1830
1831
1840FUNCTION geo_coord_inside(this, poly) RESULT(inside)
1841TYPE(geo_coord), INTENT(IN) :: this
1842TYPE(geo_coordvect), INTENT(IN) :: poly
1843LOGICAL :: inside
1844
1845INTEGER :: i, j, starti
1846
1847inside = .false.
1848IF (all(poly%ll(1,:) == poly%ll(poly%vsize,:))) THEN ! Poligono chiuso
1849 starti = 2
1850 j = 1
1851ELSE ! Poligono non chiuso
1852 starti = 1
1853 j = poly%vsize
1854ENDIF
1855DO i = starti, poly%vsize
1856 IF ((poly%ll(i,2) <= getlat(this) .AND. &
1857 getlat(this) < poly%ll(j,2)) .OR. &
1858 (poly%ll(j,2) <= getlat(this) .AND. &
1859 getlat(this) < poly%ll(i,2))) THEN
1860 IF (getlon(this) < (poly%ll(j,1) - poly%ll(i,1)) * &
1861 (getlat(this) - poly%ll(i,2)) / &
1862 (poly%ll(j,2) - poly%ll(i,2)) + poly%ll(i,1)) THEN
1864 ENDIF
1865 ENDIF
1866 j = i
1867ENDDO
1868
1869END FUNCTION geo_coord_inside
1870
1871
1872ELEMENTAL FUNCTION c_e_geo_coord(this) result (res)
1873TYPE(geo_coord),INTENT(in) :: this
1874LOGICAL :: res
1875
1876res = .not. this == geo_coord_miss
1877
1878end FUNCTION c_e_geo_coord
1879
1880
1881character(len=80) function to_char_geo_coord(this)
1882TYPE(geo_coord),INTENT(in) :: this
1883
1884to_char_geo_coord = "GEO_COORD: Lon="// &
1886 " Lat="// &
1888
1889end function to_char_geo_coord
1890
1891
1892subroutine display_geo_coord(this)
1893TYPE(geo_coord),INTENT(in) :: this
1894
1895print*,trim(to_char(this))
1896
1897end subroutine display_geo_coord
1898
1899
Export one or more geo_coordvect objects to a plain text file or to a file in ESRI/Shapefile format. Definition: geo_coord_class.F90:331 Methods for returning the value of object members. Definition: geo_coord_class.F90:293 Import one or more geo_coordvect objects from a plain text file or for a file in ESRI/Shapefile forma... Definition: geo_coord_class.F90:325 Determine whether a point lies inside a polygon or a rectangle. Definition: geo_coord_class.F90:348 Read a single geo_coord object or an array of geo_coord objects from a Fortran FORMATTED or UNFORMATT... Definition: geo_coord_class.F90:337 Represent geo_coord object in a pretty string. Definition: geo_coord_class.F90:358 Write a single geo_coord object or an array of geo_coord objects to a Fortran FORMATTED or UNFORMATTE... Definition: geo_coord_class.F90:343 Classes for handling georeferenced sparse points in geographical corodinates. Definition: geo_coord_class.F90:222 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. Definition: missing_values.f90:50 Derived type defining an isolated georeferenced point on Earth in polar geographical coordinates. Definition: geo_coord_class.F90:249 Derived type defining a one-dimensional array of georeferenced points with an associated topology (is... Definition: geo_coord_class.F90:263 |