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).
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.