libsim Versione 7.2.1

◆ geo_coordvect_exportvect()

subroutine geo_coordvect_exportvect ( type(geo_coordvect), dimension(:)  this,
character(len=*), intent(in), optional  shpfilesim,
character(len=*), intent(in), optional  shpfile,
logical, intent(in), optional  append 
)

Esporta un vettore di oggetti geo_coordvect su un file in formato testo o in formato shapefile.

Parametri
thisoggetto da esportare
[in]shpfilesimnome del file in formato testo "SIM", il parametro deve essere fornito solo se si vuole esportare su un file di quel tipo
[in]shpfilenome dello shapefile, il parametro deve essere fornito solo se si vuole esportare su un file di quel tipo
[in]appendsistema di coordinate (proiezione) dei dati, usare i parametri ::proj_geo (default) o ::proj_utm, ha senso se this, a seguito di una chiamata a ::to_geo o a ::to_utm, contiene le coordinate in entrambi i sistemi, altrimenti i dati vengono esportati automaticamente nel solo sistema disponibile
[in]appendse è presente e vale .TRUE. , ::export accoda all'eventuale file esistente anziché sovrascriverlo

Definizione alla linea 973 del file geo_coord_class.F90.

974! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
975! authors:
976! Davide Cesari <dcesari@arpa.emr.it>
977! Paolo Patruno <ppatruno@arpa.emr.it>
978
979! This program is free software; you can redistribute it and/or
980! modify it under the terms of the GNU General Public License as
981! published by the Free Software Foundation; either version 2 of
982! the License, or (at your option) any later version.
983
984! This program is distributed in the hope that it will be useful,
985! but WITHOUT ANY WARRANTY; without even the implied warranty of
986! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
987! GNU General Public License for more details.
988
989! You should have received a copy of the GNU General Public License
990! along with this program. If not, see <http://www.gnu.org/licenses/>.
991#include "config.h"
992
1001MODULE geo_coord_class
1002USE kinds
1003USE err_handling
1006!USE doubleprecision_phys_const
1008#ifdef HAVE_SHAPELIB
1009USE shapelib
1010!, ONLY: shpobject, shpreadobject, shpdestroyobject, shpwriteobject, &
1011! shpcreatesimpleobject
1012#endif
1013IMPLICIT NONE
1014
1015
1023INTEGER, PARAMETER :: fp_geo=fp_d
1024real(kind=fp_geo), PARAMETER :: fp_geo_miss=dmiss
1025
1028TYPE geo_coord
1029 PRIVATE
1030 INTEGER(kind=int_l) :: ilon, ilat
1031END TYPE geo_coord
1032
1033TYPE geo_coord_r
1034 PRIVATE
1035 REAL(kind=fp_geo) :: lon, lat
1036END TYPE geo_coord_r
1037
1038
1042TYPE geo_coordvect
1043 PRIVATE
1044 REAL(kind=fp_geo),POINTER :: ll(:,:) => null()
1045 INTEGER :: vsize, vtype
1046END TYPE geo_coordvect
1047
1049TYPE(geo_coord),PARAMETER :: geo_coord_miss=geo_coord(imiss,imiss)
1050
1051INTEGER, PARAMETER :: geo_coordvect_point = 1
1052INTEGER, PARAMETER :: geo_coordvect_arc = 3
1053INTEGER, PARAMETER :: geo_coordvect_polygon = 5
1054INTEGER, PARAMETER :: geo_coordvect_multipoint = 8
1055
1056
1060INTERFACE init
1061 MODULE PROCEDURE geo_coord_init, geo_coordvect_init
1062END INTERFACE
1063
1067INTERFACE delete
1068 MODULE PROCEDURE geo_coord_delete, geo_coordvect_delete
1069END INTERFACE
1070
1072INTERFACE getval
1073 MODULE PROCEDURE geo_coord_getval, geo_coordvect_getval
1074END INTERFACE
1075
1077INTERFACE OPERATOR (==)
1078 MODULE PROCEDURE geo_coord_eq
1079END INTERFACE
1080
1082INTERFACE OPERATOR (/=)
1083 MODULE PROCEDURE geo_coord_ne
1084END INTERFACE
1085
1086INTERFACE OPERATOR (>=)
1087 MODULE PROCEDURE geo_coord_ge
1088END INTERFACE
1089
1090INTERFACE OPERATOR (>)
1091 MODULE PROCEDURE geo_coord_gt
1092END INTERFACE
1093
1094INTERFACE OPERATOR (<=)
1095 MODULE PROCEDURE geo_coord_le
1096END INTERFACE
1097
1098INTERFACE OPERATOR (<)
1099 MODULE PROCEDURE geo_coord_lt
1100END INTERFACE
1101
1104INTERFACE import
1105 MODULE PROCEDURE geo_coordvect_import, geo_coordvect_importvect
1106END INTERFACE
1107
1110INTERFACE export
1111 MODULE PROCEDURE geo_coordvect_export, geo_coordvect_exportvect
1112END INTERFACE
1113
1116INTERFACE read_unit
1117 MODULE PROCEDURE geo_coord_read_unit, geo_coord_vect_read_unit
1118END INTERFACE
1119
1122INTERFACE write_unit
1123 MODULE PROCEDURE geo_coord_write_unit, geo_coord_vect_write_unit
1124END INTERFACE
1125
1127INTERFACE inside
1128 MODULE PROCEDURE geo_coord_inside, geo_coord_inside_rectang
1129END INTERFACE
1130
1132INTERFACE c_e
1133 MODULE PROCEDURE c_e_geo_coord
1134END INTERFACE
1135
1137INTERFACE to_char
1138 MODULE PROCEDURE to_char_geo_coord
1139END INTERFACE
1140
1142INTERFACE display
1143 MODULE PROCEDURE display_geo_coord
1144END INTERFACE
1145
1146CONTAINS
1147
1148
1149! ===================
1150! == geo_coord ==
1151! ===================
1158SUBROUTINE geo_coord_init(this, lon, lat, ilon, ilat)
1159TYPE(geo_coord) :: this
1160REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon
1161REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat
1162integer(kind=int_l), INTENT(IN), OPTIONAL :: ilon
1163integer(kind=int_l), INTENT(IN), OPTIONAL :: ilat
1164
1165real(kind=fp_geo) :: llon,llat
1166
1167CALL optio(ilon, this%ilon)
1168CALL optio(ilat, this%ilat)
1169
1170if (.not. c_e(this%ilon)) then
1171 CALL optio(lon, llon)
1172 if (c_e(llon)) this%ilon=nint(llon*1.d5)
1173end if
1174
1175if (.not. c_e(this%ilat)) then
1176 CALL optio(lat, llat)
1177 if (c_e(llat)) this%ilat=nint(llat*1.d5)
1178end if
1179
1180END SUBROUTINE geo_coord_init
1181
1183SUBROUTINE geo_coord_delete(this)
1184TYPE(geo_coord), INTENT(INOUT) :: this
1185
1186this%ilon = imiss
1187this%ilat = imiss
1188
1189END SUBROUTINE geo_coord_delete
1190
1195elemental SUBROUTINE geo_coord_getval(this, lon, lat,ilon,ilat)
1196TYPE(geo_coord),INTENT(IN) :: this
1197REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lon
1198REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lat
1199integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilon
1200integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilat
1201
1202IF (PRESENT(ilon)) ilon = getilon(this)
1203IF (PRESENT(ilat)) ilat = getilat(this)
1204
1205IF (PRESENT(lon)) lon = getlon(this)
1206IF (PRESENT(lat)) lat = getlat(this)
1207
1208END SUBROUTINE geo_coord_getval
1209
1210
1215elemental FUNCTION getilat(this)
1216TYPE(geo_coord),INTENT(IN) :: this
1217integer(kind=int_l) :: getilat
1218
1219getilat = this%ilat
1220
1221END FUNCTION getilat
1222
1227elemental FUNCTION getlat(this)
1228TYPE(geo_coord),INTENT(IN) :: this
1229real(kind=fp_geo) :: getlat
1230integer(kind=int_l) :: ilat
1231
1232ilat=getilat(this)
1233if (c_e(ilat)) then
1234 getlat = ilat*1.d-5
1235else
1236 getlat=fp_geo_miss
1237end if
1238
1239END FUNCTION getlat
1240
1245elemental FUNCTION getilon(this)
1246TYPE(geo_coord),INTENT(IN) :: this
1247integer(kind=int_l) :: getilon
1248
1249getilon = this%ilon
1250
1251END FUNCTION getilon
1252
1253
1258elemental FUNCTION getlon(this)
1259TYPE(geo_coord),INTENT(IN) :: this
1260real(kind=fp_geo) :: getlon
1261integer(kind=int_l) :: ilon
1262
1263ilon=getilon(this)
1264if (c_e(ilon)) then
1265 getlon = ilon*1.d-5
1266else
1267 getlon=fp_geo_miss
1268end if
1269
1270END FUNCTION getlon
1271
1272
1273elemental FUNCTION geo_coord_eq(this, that) RESULT(res)
1274TYPE(geo_coord),INTENT(IN) :: this, that
1275LOGICAL :: res
1276
1277res = (this%ilon == that%ilon .AND. this%ilat == that%ilat)
1278
1279END FUNCTION geo_coord_eq
1280
1281
1282elemental FUNCTION geo_coord_ne(this, that) RESULT(res)
1283TYPE(geo_coord),INTENT(IN) :: this, that
1284LOGICAL :: res
1285
1286res = .not. this == that
1287
1288END FUNCTION geo_coord_ne
1289
1292elemental FUNCTION geo_coord_gt(this, that) RESULT(res)
1293TYPE(geo_coord),INTENT(IN) :: this, that
1294LOGICAL :: res
1295
1296res = this%ilon > that%ilon
1297
1298if ( this%ilon == that%ilon ) then
1299 res = this%ilat > that%ilat
1300end if
1301
1302END FUNCTION geo_coord_gt
1303
1304
1307elemental FUNCTION geo_coord_ge(this, that) RESULT(res)
1308TYPE(geo_coord),INTENT(IN) :: this, that
1309LOGICAL :: res
1310
1311res = .not. this < that
1312
1313END FUNCTION geo_coord_ge
1314
1317elemental FUNCTION geo_coord_lt(this, that) RESULT(res)
1318TYPE(geo_coord),INTENT(IN) :: this, that
1319LOGICAL :: res
1320
1321res = this%ilon < that%ilon
1322
1323if ( this%ilon == that%ilon ) then
1324 res = this%ilat < that%ilat
1325end if
1326
1327
1328END FUNCTION geo_coord_lt
1329
1330
1333elemental FUNCTION geo_coord_le(this, that) RESULT(res)
1334TYPE(geo_coord),INTENT(IN) :: this, that
1335LOGICAL :: res
1336
1337res = .not. this > that
1338
1339END FUNCTION geo_coord_le
1340
1341
1344elemental FUNCTION geo_coord_ure(this, that) RESULT(res)
1345TYPE(geo_coord),INTENT(IN) :: this, that
1346LOGICAL :: res
1347
1348res = (this%ilon >= that%ilon .AND. this%ilat >= that%ilat)
1349
1350END FUNCTION geo_coord_ure
1351
1354elemental FUNCTION geo_coord_ur(this, that) RESULT(res)
1355TYPE(geo_coord),INTENT(IN) :: this, that
1356LOGICAL :: res
1357
1358res = (this%ilon > that%ilon .AND. this%ilat > that%ilat)
1359
1360END FUNCTION geo_coord_ur
1361
1362
1365elemental FUNCTION geo_coord_lle(this, that) RESULT(res)
1366TYPE(geo_coord),INTENT(IN) :: this, that
1367LOGICAL :: res
1368
1369res = (this%ilon <= that%ilon .AND. this%ilat <= that%ilat)
1370
1371END FUNCTION geo_coord_lle
1372
1375elemental FUNCTION geo_coord_ll(this, that) RESULT(res)
1376TYPE(geo_coord),INTENT(IN) :: this, that
1377LOGICAL :: res
1378
1379res = (this%ilon < that%ilon .AND. this%ilat < that%ilat)
1380
1381END FUNCTION geo_coord_ll
1382
1383
1389SUBROUTINE geo_coord_read_unit(this, unit)
1390TYPE(geo_coord),INTENT(out) :: this
1391INTEGER, INTENT(in) :: unit
1392
1393CALL geo_coord_vect_read_unit((/this/), unit)
1394
1395END SUBROUTINE geo_coord_read_unit
1396
1397
1403SUBROUTINE geo_coord_vect_read_unit(this, unit)
1404TYPE(geo_coord) :: this(:)
1405INTEGER, INTENT(in) :: unit
1406
1407CHARACTER(len=40) :: form
1408INTEGER :: i
1409
1410INQUIRE(unit, form=form)
1411IF (form == 'FORMATTED') THEN
1412 read(unit,*) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1413!TODO bug gfortran compiler !
1414!missing values are unredeable when formatted
1415ELSE
1416 READ(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1417ENDIF
1418
1419END SUBROUTINE geo_coord_vect_read_unit
1420
1421
1426SUBROUTINE geo_coord_write_unit(this, unit)
1427TYPE(geo_coord),INTENT(in) :: this
1428INTEGER, INTENT(in) :: unit
1429
1430CALL geo_coord_vect_write_unit((/this/), unit)
1431
1432END SUBROUTINE geo_coord_write_unit
1433
1434
1439SUBROUTINE geo_coord_vect_write_unit(this, unit)
1440TYPE(geo_coord),INTENT(in) :: this(:)
1441INTEGER, INTENT(in) :: unit
1442
1443CHARACTER(len=40) :: form
1444INTEGER :: i
1445
1446INQUIRE(unit, form=form)
1447IF (form == 'FORMATTED') THEN
1448 WRITE(unit,*) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1449!TODO bug gfortran compiler !
1450!missing values are unredeable when formatted
1451ELSE
1452 WRITE(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1453ENDIF
1454
1455END SUBROUTINE geo_coord_vect_write_unit
1456
1457
1460FUNCTION geo_coord_dist(this, that) RESULT(dist)
1462TYPE(geo_coord), INTENT (IN) :: this
1463TYPE(geo_coord), INTENT (IN) :: that
1464REAL(kind=fp_geo) :: dist
1465
1466REAL(kind=fp_geo) :: x,y
1467! Distanza approssimata, valida per piccoli angoli
1468
1469x = (getlon(this)-getlon(that))*cos(((getlat(this)+getlat(that))/2.)*degrad)
1470y = getlat(this)-getlat(that)
1471dist = sqrt(x**2 + y**2)*degrad*rearth
1472
1473END FUNCTION geo_coord_dist
1474
1475
1484FUNCTION geo_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1485TYPE(geo_coord),INTENT(IN) :: this
1486TYPE(geo_coord),INTENT(IN) :: coordmin
1487TYPE(geo_coord),INTENT(IN) :: coordmax
1488LOGICAL :: res
1489
1490res = (geo_coord_ure(this, coordmin) .AND. geo_coord_lle(this,coordmax))
1491
1492END FUNCTION geo_coord_inside_rectang
1493
1494
1495! ===================
1496! == geo_coordvect ==
1497! ===================
1506RECURSIVE SUBROUTINE geo_coordvect_init(this, lon, lat)
1507TYPE(geo_coordvect), INTENT(OUT) :: this
1508REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon(:)
1509REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat(:)
1510
1511IF (PRESENT(lon) .AND. PRESENT(lat)) THEN
1512 this%vsize = min(SIZE(lon), SIZE(lat))
1513 ALLOCATE(this%ll(this%vsize,2))
1514 this%ll(1:this%vsize,1) = lon(1:this%vsize)
1515 this%ll(1:this%vsize,2) = lat(1:this%vsize)
1516ELSE
1517 this%vsize = 0
1518 NULLIFY(this%ll)
1519ENDIF
1520this%vtype = 0 !?
1521
1522END SUBROUTINE geo_coordvect_init
1523
1524
1527SUBROUTINE geo_coordvect_delete(this)
1528TYPE(geo_coordvect), INTENT(INOUT) :: this
1529
1530IF (ASSOCIATED(this%ll)) DEALLOCATE(this%ll)
1531this%vsize = 0
1532this%vtype = 0
1533
1534END SUBROUTINE geo_coordvect_delete
1535
1536
1545SUBROUTINE geo_coordvect_getval(this, lon, lat)
1546TYPE(geo_coordvect),INTENT(IN) :: this
1547REAL(kind=fp_geo), OPTIONAL, POINTER :: lon(:)
1548REAL(kind=fp_geo), OPTIONAL, POINTER :: lat(:)
1549
1550IF (PRESENT(lon)) THEN
1551 IF (ASSOCIATED(this%ll)) THEN
1552 ALLOCATE(lon(this%vsize))
1553 lon(:) = this%ll(1:this%vsize,1)
1554 ENDIF
1555ENDIF
1556IF (PRESENT(lat)) THEN
1557 IF (ASSOCIATED(this%ll)) THEN
1558 ALLOCATE(lat(this%vsize))
1559 lat(:) = this%ll(1:this%vsize,2)
1560 ENDIF
1561ENDIF
1562
1563END SUBROUTINE geo_coordvect_getval
1564
1565
1566SUBROUTINE geo_coordvect_import(this, unitsim, &
1567#ifdef HAVE_SHAPELIB
1568 shphandle, &
1569#endif
1570 nshp)
1571TYPE(geo_coordvect), INTENT(OUT) :: this
1572INTEGER,OPTIONAL,INTENT(IN) :: unitsim
1573#ifdef HAVE_SHAPELIB
1574TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
1575#endif
1576INTEGER,OPTIONAL,INTENT(IN) :: nshp
1577
1578REAL(kind=fp_geo),ALLOCATABLE :: llon(:), llat(:)
1579REAL(kind=fp_geo) :: lv1,lv2,lv3,lv4
1580INTEGER :: i, lvsize
1581CHARACTER(len=40) :: lname
1582#ifdef HAVE_SHAPELIB
1583TYPE(shpobject) :: shpobj
1584#endif
1585
1586IF (PRESENT(unitsim)) THEN
1587 ! Leggo l'intestazione
1588 READ(unitsim,*,END=10)lvsize,lv1,lv2,lv3,lv4,lname
1589 ALLOCATE(llon(lvsize+1), llat(lvsize+1))
1590 ! Leggo il poligono
1591 READ(unitsim,*)(llon(i),llat(i), i=1,lvsize)
1592 ! Lo chiudo se necessario
1593 IF (llon(1) /= llon(lvsize) .OR. llat(1) /= llat(lvsize)) THEN
1594 lvsize = lvsize + 1
1595 llon(lvsize) = llon(1)
1596 llat(lvsize) = llat(1)
1597 ENDIF
1598 ! Lo inserisco nel mio oggetto
1599 CALL init(this, lon=llon(1:lvsize), lat=llat(1:lvsize))
1600 this%vtype = geo_coordvect_polygon ! Sempre un poligono
1601
1602 DEALLOCATE(llon, llat)
1603 RETURN
160410 CALL raise_error('nella lettura del file '//trim(to_char(unitsim)))
1605 DEALLOCATE(llon, llat) ! End of file, ritorno un oggetto non assegnato
1606#ifdef HAVE_SHAPELIB
1607ELSE IF (PRESENT(shphandle) .AND. PRESENT(nshp)) THEN
1608 ! Leggo l'oggetto shape
1609 shpobj = shpreadobject(shphandle, nshp)
1610 IF (.NOT.shpisnull(shpobj)) THEN
1611 ! Lo inserisco nel mio oggetto
1612 CALL init(this, lon=real(shpobj%padfx,kind=fp_geo), &
1613 lat=real(shpobj%padfy,kind=fp_geo))
1614 this%vtype = shpobj%nshptype
1615 CALL shpdestroyobject(shpobj)
1616 ELSE
1617 CALL init(this)
1618 ENDIF
1619#endif
1620ENDIF
1621
1622END SUBROUTINE geo_coordvect_import
1623
1624
1625SUBROUTINE geo_coordvect_export(this, unitsim, &
1626#ifdef HAVE_SHAPELIB
1627 shphandle, &
1628#endif
1629 nshp)
1630TYPE(geo_coordvect), INTENT(INOUT) :: this
1631INTEGER,OPTIONAL,INTENT(IN) :: unitsim
1632#ifdef HAVE_SHAPELIB
1633TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
1634#endif
1635INTEGER,OPTIONAL,INTENT(IN) :: nshp
1636
1637INTEGER :: i, lnshp
1638#ifdef HAVE_SHAPELIB
1639TYPE(shpobject) :: shpobj
1640#endif
1641
1642IF (PRESENT(unitsim)) THEN
1643 IF (this%vsize > 0) THEN
1644 ! Scrivo l'intestazione
1645 WRITE(unitsim,*)SIZE(this%ll,1),-1.,5000.,-0.1,1.1,'Area'
1646 ! Scrivo il poligono
1647 WRITE(unitsim,*)(this%ll(i,1:2), i=1,this%vsize)
1648 ELSE
1649 CALL raise_warning('oggetto geo_coordvect vuoto, non scrivo niente in '// &
1650 trim(to_char(unitsim)))
1651 ENDIF
1652#ifdef HAVE_SHAPELIB
1653ELSE IF (PRESENT(shphandle)) THEN
1654 IF (PRESENT(nshp)) THEN
1655 lnshp = nshp
1656 ELSE
1657 lnshp = -1 ! -1 = append
1658 ENDIF
1659 ! Creo l'oggetto shape inizializzandolo con il mio oggetto
1660 shpobj = shpcreatesimpleobject(this%vtype, this%vsize, &
1661 REAL(this%ll(1:this%vsize,1),kind=fp_d), &
1662 REAL(this%ll(1:this%vsize,2),kind=fp_d))
1663 IF (.NOT.shpisnull(shpobj)) THEN
1664 ! Lo scrivo nello shapefile
1665 i=shpwriteobject(shphandle, lnshp, shpobj)
1666 CALL shpdestroyobject(shpobj)
1667 ENDIF
1668#endif
1669ENDIF
1670
1671END SUBROUTINE geo_coordvect_export
1672
1691SUBROUTINE geo_coordvect_importvect(this, shpfilesim, shpfile)
1692TYPE(geo_coordvect),POINTER :: this(:)
1693CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim
1694CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile
1695
1696REAL(kind=fp_geo) :: inu
1697REAL(kind=fp_d) :: minb(4), maxb(4)
1698INTEGER :: i, u, ns, lvsize, shptype, dbfnf, dbfnr
1699CHARACTER(len=40) :: lname
1700#ifdef HAVE_SHAPELIB
1701TYPE(shpfileobject) :: shphandle
1702#endif
1703
1704NULLIFY(this)
1705
1706IF (PRESENT(shpfilesim)) THEN
1707 u = getunit()
1708 OPEN(u, file=shpfilesim, status='old', err=30)
1709 ns = 0 ! Conto il numero di shape contenute
1710 DO WHILE(.true.)
1711 READ(u,*,END=10,ERR=20)lvsize,inu,inu,inu,inu,lname
1712 READ(u,*,END=20,ERR=20)(inu,inu, i=1,lvsize)
1713 ns = ns + 1
1714 ENDDO
171510 CONTINUE
1716 IF (ns > 0) THEN ! Alloco e leggo il mio oggetto
1717 ALLOCATE(this(ns))
1718 rewind(u)
1719 DO i = 1, ns
1720 CALL import(this(i), unitsim=u)
1721 ENDDO
1722 ENDIF
172320 CONTINUE
1724 CLOSE(u)
1725 IF (.NOT.ASSOCIATED(this)) THEN
1726 CALL raise_warning('file '//trim(shpfilesim)//' vuoto o corrotto')
1727 ENDIF
1728 RETURN
172930 CONTINUE
1730 CALL raise_error('Impossibile aprire il file '//trim(shpfile))
1731 RETURN
1732
1733ELSE IF (PRESENT(shpfile)) THEN
1734#ifdef HAVE_SHAPELIB
1735 shphandle = shpopen(trim(shpfile), 'rb')
1736 IF (shpfileisnull(shphandle)) THEN
1737 CALL raise_error('Impossibile aprire lo shapefile '//trim(shpfile))
1738 RETURN
1739 ENDIF
1740 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr) ! Ottengo le info sul file
1741 IF (ns > 0) THEN ! Alloco e leggo il mio oggetto
1742 ALLOCATE(this(ns))
1743 this(:)%vtype = shptype
1744 DO i = 1, ns
1745 CALL import(this(i), shphandle=shphandle, nshp=i-1)
1746 ENDDO
1747 ENDIF
1748 CALL shpclose(shphandle)
1749 RETURN
1750#endif
1751ENDIF
1752
1753END SUBROUTINE geo_coordvect_importvect
1754
1755
1758SUBROUTINE geo_coordvect_exportvect(this, shpfilesim, shpfile, append)
1759TYPE(geo_coordvect) :: this(:)
1760CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim
1761CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile
1767LOGICAL,INTENT(in),OPTIONAL :: append
1768
1769REAL(kind=fp_d) :: minb(4), maxb(4)
1770INTEGER :: i, u, ns, shptype, dbfnf, dbfnr
1771LOGICAL :: lappend
1772#ifdef HAVE_SHAPELIB
1773TYPE(shpfileobject) :: shphandle
1774#endif
1775
1776IF (PRESENT(append)) THEN
1777 lappend = append
1778ELSE
1779 lappend = .false.
1780ENDIF
1781IF (PRESENT(shpfilesim)) THEN
1782 u = getunit()
1783 IF (lappend) THEN
1784 OPEN(u, file=shpfilesim, status='unknown', position='append', err=30)
1785 ELSE
1786 OPEN(u, file=shpfilesim, status='unknown', err=30)
1787 ENDIF
1788 DO i = 1, SIZE(this)
1789 CALL export(this(i), unitsim=u)
1790 ENDDO
1791 CLOSE(u)
1792 RETURN
179330 CONTINUE
1794 CALL raise_error('Impossibile aprire il file '//trim(shpfile))
1795 RETURN
1796ELSE IF (PRESENT(shpfile)) THEN
1797#ifdef HAVE_SHAPELIB
1798 IF (lappend) THEN
1799 shphandle = shpopen(trim(shpfile), 'r+b')
1800 IF (shpfileisnull(shphandle)) THEN ! shpopen funziona solo su file esistenti
1801 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
1802 ENDIF
1803 ELSE
1804 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
1805 ENDIF
1806 IF (shpfileisnull(shphandle)) THEN
1807 CALL raise_error('Impossibile aprire lo shapefile '//trim(shpfile))
1808 RETURN
1809 ENDIF
1810 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr) ! Ottengo le info sul file
1811 DO i = 1, SIZE(this)
1812 IF (i > ns .OR. lappend) THEN ! Append shape
1813 CALL export(this(i), shphandle=shphandle)
1814 ELSE ! Overwrite shape
1815 CALL export(this(i), shphandle=shphandle, nshp=i-1)
1816 ENDIF
1817 ENDDO
1818 CALL shpclose(shphandle)
1819 RETURN
1820#endif
1821ENDIF
1822
1823END SUBROUTINE geo_coordvect_exportvect
1824
1825
1834FUNCTION geo_coord_inside(this, poly) RESULT(inside)
1835TYPE(geo_coord), INTENT(IN) :: this
1836TYPE(geo_coordvect), INTENT(IN) :: poly
1837LOGICAL :: inside
1838
1839INTEGER :: i, j, starti
1840
1841inside = .false.
1842IF (all(poly%ll(1,:) == poly%ll(poly%vsize,:))) THEN ! Poligono chiuso
1843 starti = 2
1844 j = 1
1845ELSE ! Poligono non chiuso
1846 starti = 1
1847 j = poly%vsize
1848ENDIF
1849DO i = starti, poly%vsize
1850 IF ((poly%ll(i,2) <= getlat(this) .AND. &
1851 getlat(this) < poly%ll(j,2)) .OR. &
1852 (poly%ll(j,2) <= getlat(this) .AND. &
1853 getlat(this) < poly%ll(i,2))) THEN
1854 IF (getlon(this) < (poly%ll(j,1) - poly%ll(i,1)) * &
1855 (getlat(this) - poly%ll(i,2)) / &
1856 (poly%ll(j,2) - poly%ll(i,2)) + poly%ll(i,1)) THEN
1857 inside = .NOT. inside
1858 ENDIF
1859 ENDIF
1860 j = i
1861ENDDO
1862
1863END FUNCTION geo_coord_inside
1864
1865
1866ELEMENTAL FUNCTION c_e_geo_coord(this) result (res)
1867TYPE(geo_coord),INTENT(in) :: this
1868LOGICAL :: res
1869
1870res = .not. this == geo_coord_miss
1871
1872end FUNCTION c_e_geo_coord
1873
1874
1875character(len=80) function to_char_geo_coord(this)
1876TYPE(geo_coord),INTENT(in) :: this
1877
1878to_char_geo_coord = "GEO_COORD: Lon="// &
1879 trim(to_char(getlon(this),miss="Missing lon",form="(f11.5)"))//&
1880 " Lat="// &
1881 trim(to_char(getlat(this),miss="Missing lat",form="(f11.5)"))
1882
1883end function to_char_geo_coord
1884
1885
1886subroutine display_geo_coord(this)
1887TYPE(geo_coord),INTENT(in) :: this
1888
1889print*,trim(to_char(this))
1890
1891end subroutine display_geo_coord
1892
1893
1894END 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.