libsim Versione 7.2.6

◆ 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!> Classes for handling georeferenced sparse points in geographical
993!! corodinates. This module defines two classes for managing
994!! georeferenced points on the Earth in geographical polar
995!! coordinates. It allows importing and exporting blocks of points
996!! from/to plain text files and \a ESRI \a shapefile's. Both classes
997!! have \a PRIVATE members, so that they cannot be manipulated
998!! directly, but only through the proper methods.
999!!
1000!! \ingroup base
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
1016!> REAL Kind for geographical coordinates.
1017!! This constant has to be used when defining or converting values to be
1018!! used as geographical coordinates, e.g.:
1019!! \code
1020!! REAL(kind=fp_geo) :: x, y
1021!! coordx = REAL(mylon, kind=fp_geo)
1022!! \endcode
1023INTEGER, PARAMETER :: fp_geo=fp_d
1024real(kind=fp_geo), PARAMETER :: fp_geo_miss=dmiss
1025
1026!> Derived type defining an isolated georeferenced point on Earth in
1027!! polar geographical coordinates
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
1039!> Derived type defining a one-dimensional array of georeferenced points
1040!! with an associated topology (isolated point, arc, polygon, group of
1041!! points)
1042TYPE geo_coordvect
1043 PRIVATE
1044 REAL(kind=fp_geo),POINTER :: ll(:,:) => null()
1045 INTEGER :: vsize, vtype
1046END TYPE geo_coordvect
1047
1048!> Missing value for geo_coord
1049TYPE(geo_coord),PARAMETER :: geo_coord_miss=geo_coord(imiss,imiss)
1050
1051INTEGER, PARAMETER :: geo_coordvect_point = 1 !< Topology for geo_coordvect (from shapelib): isolated point
1052INTEGER, PARAMETER :: geo_coordvect_arc = 3 !< Topology for geo_coordvect (from shapelib): arc (multiple arcs unsupported)
1053INTEGER, PARAMETER :: geo_coordvect_polygon = 5 !< Topology for geo_coordvect (from shapelib): polygon (necessarily closed, multiple polygons unsupported)
1054INTEGER, PARAMETER :: geo_coordvect_multipoint = 8 !< Topology for geo_coordvect (from shapelib): group of points
1055
1056
1057!> Constructors for the two classes.
1058!! They have to be called for every object of these types which is
1059!! going to be used.
1060INTERFACE init
1061 MODULE PROCEDURE geo_coord_init, geo_coordvect_init
1062END INTERFACE
1063
1064!> Detructors for the two classes.
1065!! They clean up all the information associated with the corresponding
1066!! objects.
1067INTERFACE delete
1068 MODULE PROCEDURE geo_coord_delete, geo_coordvect_delete
1069END INTERFACE
1070
1071!> Methods for returning the value of object members.
1072INTERFACE getval
1073 MODULE PROCEDURE geo_coord_getval, geo_coordvect_getval
1074END INTERFACE
1075
1076!> Logical equality operator.
1077INTERFACE OPERATOR (==)
1078 MODULE PROCEDURE geo_coord_eq
1079END INTERFACE
1080
1081!> Logical inequality operator.
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
1102!> Import one or more \a geo_coordvect objects from a plain text file
1103!! or for a file in ESRI/Shapefile format.
1104INTERFACE import
1105 MODULE PROCEDURE geo_coordvect_import, geo_coordvect_importvect
1106END INTERFACE
1107
1108!> Export one or more \a geo_coordvect objects to a plain text file
1109!! or to a file in ESRI/Shapefile format.
1110INTERFACE export
1111 MODULE PROCEDURE geo_coordvect_export, geo_coordvect_exportvect
1112END INTERFACE
1113
1114!> Read a single \a geo_coord object or an array of \a geo_coord objects
1115!! from a Fortran \c FORMATTED or \c UNFORMATTED file.
1116INTERFACE read_unit
1117 MODULE PROCEDURE geo_coord_read_unit, geo_coord_vect_read_unit
1118END INTERFACE
1119
1120!> Write a single \a geo_coord object or an array of \a geo_coord objects
1121!! to a Fortran \c FORMATTED or \c UNFORMATTED file.
1122INTERFACE write_unit
1123 MODULE PROCEDURE geo_coord_write_unit, geo_coord_vect_write_unit
1124END INTERFACE
1125
1126!> Determine whether a point lies inside a polygon or a rectangle.
1127INTERFACE inside
1128 MODULE PROCEDURE geo_coord_inside, geo_coord_inside_rectang
1129END INTERFACE
1130
1131!> Missing check
1132INTERFACE c_e
1133 MODULE PROCEDURE c_e_geo_coord
1134END INTERFACE
1135
1136!>Represent geo_coord object in a pretty string
1137INTERFACE to_char
1138 MODULE PROCEDURE to_char_geo_coord
1139END INTERFACE
1140
1141!>Print object
1142INTERFACE display
1143 MODULE PROCEDURE display_geo_coord
1144END INTERFACE
1145
1146CONTAINS
1147
1148
1149! ===================
1150! == geo_coord ==
1151! ===================
1152!> Costruisce un oggetto \a geo_coord con i parametri opzionali forniti.
1153!! Se sono presenti \a lon e \a lat, inizializza le coordinate geografiche
1154!! ignorando \a utme e \a utmn, mentre se sono specificati \a utme e \a utmn
1155!! succede il contrario; non è possibile specificare le coordinate in entrambi
1156!! i sistemi, usare eventualmente \a to_geo. Se non viene passato nessun parametro
1157!! opzionale l'oggetto è inizializzato a valore mancante.
1158SUBROUTINE geo_coord_init(this, lon, lat, ilon, ilat)
1159TYPE(geo_coord) :: this !< oggetto da inizializzare
1160REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon !< longitudine geografica
1161REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat !< latitudine geografica
1162integer(kind=int_l), INTENT(IN), OPTIONAL :: ilon !< integer longitudine geografica (nint(lon*1.e5)
1163integer(kind=int_l), INTENT(IN), OPTIONAL :: ilat !< integer latitudine geografica (nint(lat*1.e5)
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
1182!> Distrugge l'oggetto in maniera pulita, assegnandogli un valore mancante.
1183SUBROUTINE geo_coord_delete(this)
1184TYPE(geo_coord), INTENT(INOUT) :: this !< oggetto da distruggre
1185
1186this%ilon = imiss
1187this%ilat = imiss
1188
1189END SUBROUTINE geo_coord_delete
1190
1191!> Restituisce il valore di uno o più componenti di un oggetto \a geo_coord.
1192!! Qualsiasi combinazione dei parametri opzionali è consentita; se
1193!! il tipo di coordinata richiesta non è stato inizializzato né calcolato,
1194!! restituisce il corrispondente valore mancante.
1195elemental SUBROUTINE geo_coord_getval(this, lon, lat,ilon,ilat)
1196TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire i componenti
1197REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lon !< longitudine geografica
1198REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lat !< latitudine geografica
1199integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilon !< integer longitudine geografica (nint(lon*1.e5)
1200integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilat !< integer latitudine geografica (nint(lat*1.e5)
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
1211!> Restituisce la latitudine di uno o più componenti di un oggetto \a geo_coord.
1212!! Se la latitudine non è stata inizializzata né calcolata
1213!! restituisce il corrispondente valore mancante.
1214!! Nata per permettere operazioni vettorizzate
1215elemental FUNCTION getilat(this)
1216TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire latitudine
1217integer(kind=int_l) :: getilat !< latitudine geografica
1218
1219getilat = this%ilat
1220
1221END FUNCTION getilat
1222
1223!> Restituisce la latitudine di uno o più componenti di un oggetto \a geo_coord.
1224!! Se la latitudine non è stata inizializzata né calcolata
1225!! restituisce il corrispondente valore mancante.
1226!! Nata per permettere operazioni vettorizzate
1227elemental FUNCTION getlat(this)
1228TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire latitudine
1229real(kind=fp_geo) :: getlat !< latitudine geografica
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
1241!> Restituisce la longitudine di uno o più componenti di un oggetto \a geo_coord.
1242!! Se la latitudine non è stata inizializzata né calcolata
1243!! restituisce il corrispondente valore mancante.
1244!! Nata per permettere operazioni vettorizzate
1245elemental FUNCTION getilon(this)
1246TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire latitudine
1247integer(kind=int_l) :: getilon !< longitudine geografica
1248
1249getilon = this%ilon
1250
1251END FUNCTION getilon
1252
1253
1254!> Restituisce la longitudine di uno o più componenti di un oggetto \a geo_coord.
1255!! Se la latitudine non è stata inizializzata né calcolata
1256!! restituisce il corrispondente valore mancante.
1257!! Nata per permettere operazioni vettorizzate
1258elemental FUNCTION getlon(this)
1259TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire latitudine
1260real(kind=fp_geo) :: getlon !< longitudine geografica
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
1290!> Logical great operator. Returns true if the first point
1291!! lies to the west of the second or if lon is equal north
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
1305!> Logical great-equal operator. Returns true if the first point
1306!! lies to the west of the second or if lon is equal north
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
1315!> Logical less operator. Returns true if the first point
1316!! lies to the west of the second or if lon is equal south
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
1331!> Logical less-equal operator. Returns true if the first point
1332!! lies to the west of the second or if lon is equal south
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
1342!> Logical greater-equal operator. Returns true if the first point
1343!! lies to the northeast of the second
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
1352!> Logical greater operator. Returns true if the first point
1353!! lies to the northeast of the second
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
1363!> Logical less-equal operator. Returns true if the first point
1364!! lies to the southwest of the second
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
1373!> Logical less operator. Returns true if the first point
1374!! lies to the southwest of the second
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
1384!> Legge da un'unità di file il contenuto dell'oggetto \a this.
1385!! Il record da leggere deve essere stato scritto con la ::write_unit
1386!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
1387!! del vettore devono essere accordate. Il metodo controlla se il file è
1388!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1389SUBROUTINE geo_coord_read_unit(this, unit)
1390TYPE(geo_coord),INTENT(out) :: this !< oggetto da leggere
1391INTEGER, INTENT(in) :: unit !< unità da cui leggere
1392
1393CALL geo_coord_vect_read_unit((/this/), unit)
1394
1395END SUBROUTINE geo_coord_read_unit
1396
1397
1398!> Legge da un'unità di file il contenuto dell'oggetto \a this.
1399!! Il record da leggere deve essere stato scritto con la ::write_unit
1400!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
1401!! del vettore devono essere accordate. Il metodo controlla se il file è
1402!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1403SUBROUTINE geo_coord_vect_read_unit(this, unit)
1404TYPE(geo_coord) :: this(:) !< oggetto da leggere
1405INTEGER, INTENT(in) :: unit !< unità da cui leggere
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
1422!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
1423!! Il record scritto potrà successivamente essere letto con la ::read_unit.
1424!! Il metodo controlla se il file è
1425!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1426SUBROUTINE geo_coord_write_unit(this, unit)
1427TYPE(geo_coord),INTENT(in) :: this !< oggetto da scrivere
1428INTEGER, INTENT(in) :: unit !< unità su cui scrivere
1429
1430CALL geo_coord_vect_write_unit((/this/), unit)
1431
1432END SUBROUTINE geo_coord_write_unit
1433
1434
1435!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
1436!! Il record scritto potrà successivamente essere letto con la ::read_unit.
1437!! Il metodo controlla se il file è
1438!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1439SUBROUTINE geo_coord_vect_write_unit(this, unit)
1440TYPE(geo_coord),INTENT(in) :: this(:) !< oggetto da scrivere
1441INTEGER, INTENT(in) :: unit !< unità su cui scrivere
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
1458!> Restituisce la distanza in m tra 2 oggetti geo_coord.
1459!! La distanza è calcolata approssimativamente ed è valida per piccoli angoli.
1460FUNCTION geo_coord_dist(this, that) RESULT(dist)
1462TYPE(geo_coord), INTENT (IN) :: this !< primo punto
1463TYPE(geo_coord), INTENT (IN) :: that !< secondo punto
1464REAL(kind=fp_geo) :: dist !< distanza in metri
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
1476!> Determina se il punto indicato da \a this è contenuto in un rettangolo.
1477!! Il rettangolo è orientato parallelamente agli assi del sistema,
1478!! i suoi vertici sud-ovest e nord-est sono specificati da altri due punti.
1479!! La funzione restituisce \c .TRUE. anche se il punto si trova sulla frontiera
1480!! del rettangolo.
1481!! Tutti gli oggetti devono essere già stati
1482!! convertiti ad un sistema di coordinate comune, altrimenti viene
1483!! restituito \c .FALSE. .
1484FUNCTION geo_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1485TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui determinare la posizione
1486TYPE(geo_coord),INTENT(IN) :: coordmin !< vertice sud-ovest del rettangolo
1487TYPE(geo_coord),INTENT(IN) :: coordmax !< vertice nord-est del rettangolo
1488LOGICAL :: res !< \c .TRUE. se \a this è dentro il rettangolo o sul bordo e \c .FALSE. se è fuori
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! ===================
1498!> Costruisce un oggetto \a geo_coordvect con i parametri opzionali forniti.
1499!! Se sono presenti \a lon e \a lat, inizializza le coordinate geografiche
1500!! ignorando \a utme e \a utmn, mentre se sono specificati \a utme e \a utmn
1501!! succede il contrario; non è possibile specificare le coordinate in entrambi
1502!! i sistemi, usare eventualmente \a to_geo. Se non viene passato nessun parametro
1503!! opzionale l'oggetto è inizializzato a valore mancante.
1504!! Il numero di punti dell'oggetto finale sarà uguale all'estensione
1505!! del più breve vettore della coppia fornita.
1506RECURSIVE SUBROUTINE geo_coordvect_init(this, lon, lat)
1507TYPE(geo_coordvect), INTENT(OUT) :: this !< oggetto da inizializzare
1508REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon(:) !< longitudine geografica
1509REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat(:) !< latitudine geografica
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
1525!> Distrugge l'oggetto in maniera pulita, liberando l'eventuale spazio
1526!! dinamicamente allocato.
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
1537!> Restituisce il valore di uno o più componenti di un oggetto \a geo_coordvect.
1538!! Qualsiasi combinazione dei parametri opzionali è consentita; se
1539!! il tipo di coordinata richiesta non è stato inizializzato né calcolato,
1540!! restituisce il corrispondente valore mancante.
1541!! Se forniti, i parametri \a lon, \a lat, \a utme, \a utmn devono essere
1542!! dichiarati come puntatori che vengono
1543!! allocati dalla \a getval stessa e che devono poi essere
1544!! deallocati esplicitamente dal programma chiamante.
1545SUBROUTINE geo_coordvect_getval(this, lon, lat)
1546TYPE(geo_coordvect),INTENT(IN) :: this !< oggetto di cui restituire i componenti
1547REAL(kind=fp_geo), OPTIONAL, POINTER :: lon(:) !< longitudine geografica
1548REAL(kind=fp_geo), OPTIONAL, POINTER :: lat(:) !< latitudine geografica
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
1673!> Importa un vettore di oggetti \a geo_coordvect da un file in
1674!! formato testo o in formato \a shapefile.
1675!! Il parametro \a this è un puntatore che sarà allocato a cura del metodo stesso e
1676!! dovrà invece essere deallocato da parte del programma chiamante
1677!! dopo aver chiamato il metodo \a delete per ogni suo elemento.
1678!! In caso di errore nella fase iniziale di importazione,
1679!! \a this non verrà associato, e quindi è opportuno testare
1680!! \code
1681!! IF (ASSOCIATED(my_coord_vect)) THEN...
1682!! \endcode
1683!! nel programma chiamante
1684!! per intrappolare eventuale condizioni di errore (tipicamente file non
1685!! trovato o in un formato non compatibile).
1686!! Entrambi i formati di ingresso non contengono informazioni sul tipo
1687!! di coordinate dei dati
1688!! (per il formato shapefile è possibile solo con delle estensioni non standard),
1689!! per cui questa informazione, se desiderata, deve essere fornita dal
1690!! programma chiamante.
1691SUBROUTINE geo_coordvect_importvect(this, shpfilesim, shpfile)
1692TYPE(geo_coordvect),POINTER :: this(:) !< puntatore all'oggetto su cui importare i dati, viene allocato dalla \a import stessa
1693CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim !< nome del file in formato testo "SIM", il parametro deve essere fornito solo se si vuole importare da un file di quel tipo
1694CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile !< nome delllo shapefile, il parametro deve essere fornito solo se si vuole importare da un file di quel tipo
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
1756!> Esporta un vettore di oggetti \a geo_coordvect su un file in
1757!! formato testo o in formato \a shapefile.
1758SUBROUTINE geo_coordvect_exportvect(this, shpfilesim, shpfile, append)
1759TYPE(geo_coordvect) :: this(:) !< oggetto da esportare
1760CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim !< nome del file in formato testo "SIM", il parametro deve essere fornito solo se si vuole esportare su un file di quel tipo
1761CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile !< nome dello shapefile, il parametro deve essere fornito solo se si vuole esportare su un file di quel tipo
1762!> sistema di coordinate (proiezione) dei dati,
1763!! usare i parametri \a ::proj_geo (default) o \a ::proj_utm,
1764!! ha senso se \a this, a seguito di una chiamata a \a ::to_geo o a \a ::to_utm,
1765!! contiene le coordinate in entrambi i sistemi, altrimenti i dati vengono
1766!! esportati automaticamente nel solo sistema disponibile
1767LOGICAL,INTENT(in),OPTIONAL :: append !< se è presente e vale \c .TRUE. , ::export accoda all'eventuale file esistente anziché sovrascriverlo
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
1826!> Determina se il punto indicato da \a this si trova
1827!! dentro o fuori dal poligono descritto dall'oggetto \a poly.
1828!! Funziona anche se la topologia di \a poly non è poligonale,
1829!! forzandone la chiusura; usa un algoritmo di ricerca del numero di
1830!! intersezioni, come indicato in
1831!! comp.graphics.algorithms FAQ (http://www.faqs.org/faqs/graphics/algorithms-faq/)
1832!! o in
1833!! http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
1834FUNCTION geo_coord_inside(this, poly) RESULT(inside)
1835TYPE(geo_coord), INTENT(IN) :: this !< oggetto di cui determinare la posizione
1836TYPE(geo_coordvect), INTENT(IN) :: poly !< poligono contenitore
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.