libsim Versione 7.2.6

◆ geo_coord_inside()

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

Determina se il punto indicato da this si trova dentro o fuori dal poligono descritto dall'oggetto poly.

Funziona anche se la topologia di poly non è poligonale, forzandone la chiusura; usa un algoritmo di ricerca del numero di intersezioni, come indicato in comp.graphics.algorithms FAQ (http://www.faqs.org/faqs/graphics/algorithms-faq/) o in http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

Parametri
[in]thisoggetto di cui determinare la posizione
[in]polypoligono contenitore

Definizione alla linea 1049 del file geo_coord_class.F90.

1050! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1051! authors:
1052! Davide Cesari <dcesari@arpa.emr.it>
1053! Paolo Patruno <ppatruno@arpa.emr.it>
1054
1055! This program is free software; you can redistribute it and/or
1056! modify it under the terms of the GNU General Public License as
1057! published by the Free Software Foundation; either version 2 of
1058! the License, or (at your option) any later version.
1059
1060! This program is distributed in the hope that it will be useful,
1061! but WITHOUT ANY WARRANTY; without even the implied warranty of
1062! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1063! GNU General Public License for more details.
1064
1065! You should have received a copy of the GNU General Public License
1066! along with this program. If not, see <http://www.gnu.org/licenses/>.
1067#include "config.h"
1068!> Classes for handling georeferenced sparse points in geographical
1069!! corodinates. This module defines two classes for managing
1070!! georeferenced points on the Earth in geographical polar
1071!! coordinates. It allows importing and exporting blocks of points
1072!! from/to plain text files and \a ESRI \a shapefile's. Both classes
1073!! have \a PRIVATE members, so that they cannot be manipulated
1074!! directly, but only through the proper methods.
1075!!
1076!! \ingroup base
1077MODULE geo_coord_class
1078USE kinds
1079USE err_handling
1082!USE doubleprecision_phys_const
1084#ifdef HAVE_SHAPELIB
1085USE shapelib
1086!, ONLY: shpobject, shpreadobject, shpdestroyobject, shpwriteobject, &
1087! shpcreatesimpleobject
1088#endif
1089IMPLICIT NONE
1090
1091
1092!> REAL Kind for geographical coordinates.
1093!! This constant has to be used when defining or converting values to be
1094!! used as geographical coordinates, e.g.:
1095!! \code
1096!! REAL(kind=fp_geo) :: x, y
1097!! coordx = REAL(mylon, kind=fp_geo)
1098!! \endcode
1099INTEGER, PARAMETER :: fp_geo=fp_d
1100real(kind=fp_geo), PARAMETER :: fp_geo_miss=dmiss
1101
1102!> Derived type defining an isolated georeferenced point on Earth in
1103!! polar geographical coordinates
1104TYPE geo_coord
1105 PRIVATE
1106 INTEGER(kind=int_l) :: ilon, ilat
1107END TYPE geo_coord
1108
1109TYPE geo_coord_r
1110 PRIVATE
1111 REAL(kind=fp_geo) :: lon, lat
1112END TYPE geo_coord_r
1113
1114
1115!> Derived type defining a one-dimensional array of georeferenced points
1116!! with an associated topology (isolated point, arc, polygon, group of
1117!! points)
1118TYPE geo_coordvect
1119 PRIVATE
1120 REAL(kind=fp_geo),POINTER :: ll(:,:) => null()
1121 INTEGER :: vsize, vtype
1122END TYPE geo_coordvect
1123
1124!> Missing value for geo_coord
1125TYPE(geo_coord),PARAMETER :: geo_coord_miss=geo_coord(imiss,imiss)
1126
1127INTEGER, PARAMETER :: geo_coordvect_point = 1 !< Topology for geo_coordvect (from shapelib): isolated point
1128INTEGER, PARAMETER :: geo_coordvect_arc = 3 !< Topology for geo_coordvect (from shapelib): arc (multiple arcs unsupported)
1129INTEGER, PARAMETER :: geo_coordvect_polygon = 5 !< Topology for geo_coordvect (from shapelib): polygon (necessarily closed, multiple polygons unsupported)
1130INTEGER, PARAMETER :: geo_coordvect_multipoint = 8 !< Topology for geo_coordvect (from shapelib): group of points
1131
1132
1133!> Constructors for the two classes.
1134!! They have to be called for every object of these types which is
1135!! going to be used.
1136INTERFACE init
1137 MODULE PROCEDURE geo_coord_init, geo_coordvect_init
1138END INTERFACE
1139
1140!> Detructors for the two classes.
1141!! They clean up all the information associated with the corresponding
1142!! objects.
1143INTERFACE delete
1144 MODULE PROCEDURE geo_coord_delete, geo_coordvect_delete
1145END INTERFACE
1146
1147!> Methods for returning the value of object members.
1148INTERFACE getval
1149 MODULE PROCEDURE geo_coord_getval, geo_coordvect_getval
1150END INTERFACE
1151
1152!> Logical equality operator.
1153INTERFACE OPERATOR (==)
1154 MODULE PROCEDURE geo_coord_eq
1155END INTERFACE
1156
1157!> Logical inequality operator.
1158INTERFACE OPERATOR (/=)
1159 MODULE PROCEDURE geo_coord_ne
1160END INTERFACE
1161
1162INTERFACE OPERATOR (>=)
1163 MODULE PROCEDURE geo_coord_ge
1164END INTERFACE
1165
1166INTERFACE OPERATOR (>)
1167 MODULE PROCEDURE geo_coord_gt
1168END INTERFACE
1169
1170INTERFACE OPERATOR (<=)
1171 MODULE PROCEDURE geo_coord_le
1172END INTERFACE
1173
1174INTERFACE OPERATOR (<)
1175 MODULE PROCEDURE geo_coord_lt
1176END INTERFACE
1177
1178!> Import one or more \a geo_coordvect objects from a plain text file
1179!! or for a file in ESRI/Shapefile format.
1180INTERFACE import
1181 MODULE PROCEDURE geo_coordvect_import, geo_coordvect_importvect
1182END INTERFACE
1183
1184!> Export one or more \a geo_coordvect objects to a plain text file
1185!! or to a file in ESRI/Shapefile format.
1186INTERFACE export
1187 MODULE PROCEDURE geo_coordvect_export, geo_coordvect_exportvect
1188END INTERFACE
1189
1190!> Read a single \a geo_coord object or an array of \a geo_coord objects
1191!! from a Fortran \c FORMATTED or \c UNFORMATTED file.
1192INTERFACE read_unit
1193 MODULE PROCEDURE geo_coord_read_unit, geo_coord_vect_read_unit
1194END INTERFACE
1195
1196!> Write a single \a geo_coord object or an array of \a geo_coord objects
1197!! to a Fortran \c FORMATTED or \c UNFORMATTED file.
1198INTERFACE write_unit
1199 MODULE PROCEDURE geo_coord_write_unit, geo_coord_vect_write_unit
1200END INTERFACE
1201
1202!> Determine whether a point lies inside a polygon or a rectangle.
1203INTERFACE inside
1204 MODULE PROCEDURE geo_coord_inside, geo_coord_inside_rectang
1205END INTERFACE
1206
1207!> Missing check
1208INTERFACE c_e
1209 MODULE PROCEDURE c_e_geo_coord
1210END INTERFACE
1211
1212!>Represent geo_coord object in a pretty string
1213INTERFACE to_char
1214 MODULE PROCEDURE to_char_geo_coord
1215END INTERFACE
1216
1217!>Print object
1218INTERFACE display
1219 MODULE PROCEDURE display_geo_coord
1220END INTERFACE
1221
1222CONTAINS
1223
1224
1225! ===================
1226! == geo_coord ==
1227! ===================
1228!> Costruisce un oggetto \a geo_coord con i parametri opzionali forniti.
1229!! Se sono presenti \a lon e \a lat, inizializza le coordinate geografiche
1230!! ignorando \a utme e \a utmn, mentre se sono specificati \a utme e \a utmn
1231!! succede il contrario; non è possibile specificare le coordinate in entrambi
1232!! i sistemi, usare eventualmente \a to_geo. Se non viene passato nessun parametro
1233!! opzionale l'oggetto è inizializzato a valore mancante.
1234SUBROUTINE geo_coord_init(this, lon, lat, ilon, ilat)
1235TYPE(geo_coord) :: this !< oggetto da inizializzare
1236REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon !< longitudine geografica
1237REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat !< latitudine geografica
1238integer(kind=int_l), INTENT(IN), OPTIONAL :: ilon !< integer longitudine geografica (nint(lon*1.e5)
1239integer(kind=int_l), INTENT(IN), OPTIONAL :: ilat !< integer latitudine geografica (nint(lat*1.e5)
1240
1241real(kind=fp_geo) :: llon,llat
1242
1243CALL optio(ilon, this%ilon)
1244CALL optio(ilat, this%ilat)
1245
1246if (.not. c_e(this%ilon)) then
1247 CALL optio(lon, llon)
1248 if (c_e(llon)) this%ilon=nint(llon*1.d5)
1249end if
1250
1251if (.not. c_e(this%ilat)) then
1252 CALL optio(lat, llat)
1253 if (c_e(llat)) this%ilat=nint(llat*1.d5)
1254end if
1255
1256END SUBROUTINE geo_coord_init
1257
1258!> Distrugge l'oggetto in maniera pulita, assegnandogli un valore mancante.
1259SUBROUTINE geo_coord_delete(this)
1260TYPE(geo_coord), INTENT(INOUT) :: this !< oggetto da distruggre
1261
1262this%ilon = imiss
1263this%ilat = imiss
1264
1265END SUBROUTINE geo_coord_delete
1266
1267!> Restituisce il valore di uno o più componenti di un oggetto \a geo_coord.
1268!! Qualsiasi combinazione dei parametri opzionali è consentita; se
1269!! il tipo di coordinata richiesta non è stato inizializzato né calcolato,
1270!! restituisce il corrispondente valore mancante.
1271elemental SUBROUTINE geo_coord_getval(this, lon, lat,ilon,ilat)
1272TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire i componenti
1273REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lon !< longitudine geografica
1274REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lat !< latitudine geografica
1275integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilon !< integer longitudine geografica (nint(lon*1.e5)
1276integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilat !< integer latitudine geografica (nint(lat*1.e5)
1277
1278IF (PRESENT(ilon)) ilon = getilon(this)
1279IF (PRESENT(ilat)) ilat = getilat(this)
1280
1281IF (PRESENT(lon)) lon = getlon(this)
1282IF (PRESENT(lat)) lat = getlat(this)
1283
1284END SUBROUTINE geo_coord_getval
1285
1286
1287!> Restituisce la latitudine di uno o più componenti di un oggetto \a geo_coord.
1288!! Se la latitudine non è stata inizializzata né calcolata
1289!! restituisce il corrispondente valore mancante.
1290!! Nata per permettere operazioni vettorizzate
1291elemental FUNCTION getilat(this)
1292TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire latitudine
1293integer(kind=int_l) :: getilat !< latitudine geografica
1294
1295getilat = this%ilat
1296
1297END FUNCTION getilat
1298
1299!> Restituisce la latitudine di uno o più componenti di un oggetto \a geo_coord.
1300!! Se la latitudine non è stata inizializzata né calcolata
1301!! restituisce il corrispondente valore mancante.
1302!! Nata per permettere operazioni vettorizzate
1303elemental FUNCTION getlat(this)
1304TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire latitudine
1305real(kind=fp_geo) :: getlat !< latitudine geografica
1306integer(kind=int_l) :: ilat
1307
1308ilat=getilat(this)
1309if (c_e(ilat)) then
1310 getlat = ilat*1.d-5
1311else
1312 getlat=fp_geo_miss
1313end if
1314
1315END FUNCTION getlat
1316
1317!> Restituisce la longitudine di uno o più componenti di un oggetto \a geo_coord.
1318!! Se la latitudine non è stata inizializzata né calcolata
1319!! restituisce il corrispondente valore mancante.
1320!! Nata per permettere operazioni vettorizzate
1321elemental FUNCTION getilon(this)
1322TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire latitudine
1323integer(kind=int_l) :: getilon !< longitudine geografica
1324
1325getilon = this%ilon
1326
1327END FUNCTION getilon
1328
1329
1330!> Restituisce la longitudine di uno o più componenti di un oggetto \a geo_coord.
1331!! Se la latitudine non è stata inizializzata né calcolata
1332!! restituisce il corrispondente valore mancante.
1333!! Nata per permettere operazioni vettorizzate
1334elemental FUNCTION getlon(this)
1335TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire latitudine
1336real(kind=fp_geo) :: getlon !< longitudine geografica
1337integer(kind=int_l) :: ilon
1338
1339ilon=getilon(this)
1340if (c_e(ilon)) then
1341 getlon = ilon*1.d-5
1342else
1343 getlon=fp_geo_miss
1344end if
1345
1346END FUNCTION getlon
1347
1348
1349elemental FUNCTION geo_coord_eq(this, that) RESULT(res)
1350TYPE(geo_coord),INTENT(IN) :: this, that
1351LOGICAL :: res
1352
1353res = (this%ilon == that%ilon .AND. this%ilat == that%ilat)
1354
1355END FUNCTION geo_coord_eq
1356
1357
1358elemental FUNCTION geo_coord_ne(this, that) RESULT(res)
1359TYPE(geo_coord),INTENT(IN) :: this, that
1360LOGICAL :: res
1361
1362res = .not. this == that
1363
1364END FUNCTION geo_coord_ne
1365
1366!> Logical great operator. Returns true if the first point
1367!! lies to the west of the second or if lon is equal north
1368elemental FUNCTION geo_coord_gt(this, that) RESULT(res)
1369TYPE(geo_coord),INTENT(IN) :: this, that
1370LOGICAL :: res
1371
1372res = this%ilon > that%ilon
1373
1374if ( this%ilon == that%ilon ) then
1375 res = this%ilat > that%ilat
1376end if
1377
1378END FUNCTION geo_coord_gt
1379
1380
1381!> Logical great-equal operator. Returns true if the first point
1382!! lies to the west of the second or if lon is equal north
1383elemental FUNCTION geo_coord_ge(this, that) RESULT(res)
1384TYPE(geo_coord),INTENT(IN) :: this, that
1385LOGICAL :: res
1386
1387res = .not. this < that
1388
1389END FUNCTION geo_coord_ge
1390
1391!> Logical less operator. Returns true if the first point
1392!! lies to the west of the second or if lon is equal south
1393elemental FUNCTION geo_coord_lt(this, that) RESULT(res)
1394TYPE(geo_coord),INTENT(IN) :: this, that
1395LOGICAL :: res
1396
1397res = this%ilon < that%ilon
1398
1399if ( this%ilon == that%ilon ) then
1400 res = this%ilat < that%ilat
1401end if
1402
1403
1404END FUNCTION geo_coord_lt
1405
1406
1407!> Logical less-equal operator. Returns true if the first point
1408!! lies to the west of the second or if lon is equal south
1409elemental FUNCTION geo_coord_le(this, that) RESULT(res)
1410TYPE(geo_coord),INTENT(IN) :: this, that
1411LOGICAL :: res
1412
1413res = .not. this > that
1414
1415END FUNCTION geo_coord_le
1416
1417
1418!> Logical greater-equal operator. Returns true if the first point
1419!! lies to the northeast of the second
1420elemental FUNCTION geo_coord_ure(this, that) RESULT(res)
1421TYPE(geo_coord),INTENT(IN) :: this, that
1422LOGICAL :: res
1423
1424res = (this%ilon >= that%ilon .AND. this%ilat >= that%ilat)
1425
1426END FUNCTION geo_coord_ure
1427
1428!> Logical greater operator. Returns true if the first point
1429!! lies to the northeast of the second
1430elemental FUNCTION geo_coord_ur(this, that) RESULT(res)
1431TYPE(geo_coord),INTENT(IN) :: this, that
1432LOGICAL :: res
1433
1434res = (this%ilon > that%ilon .AND. this%ilat > that%ilat)
1435
1436END FUNCTION geo_coord_ur
1437
1438
1439!> Logical less-equal operator. Returns true if the first point
1440!! lies to the southwest of the second
1441elemental FUNCTION geo_coord_lle(this, that) RESULT(res)
1442TYPE(geo_coord),INTENT(IN) :: this, that
1443LOGICAL :: res
1444
1445res = (this%ilon <= that%ilon .AND. this%ilat <= that%ilat)
1446
1447END FUNCTION geo_coord_lle
1448
1449!> Logical less operator. Returns true if the first point
1450!! lies to the southwest of the second
1451elemental FUNCTION geo_coord_ll(this, that) RESULT(res)
1452TYPE(geo_coord),INTENT(IN) :: this, that
1453LOGICAL :: res
1454
1455res = (this%ilon < that%ilon .AND. this%ilat < that%ilat)
1456
1457END FUNCTION geo_coord_ll
1458
1459
1460!> Legge da un'unità di file il contenuto dell'oggetto \a this.
1461!! Il record da leggere deve essere stato scritto con la ::write_unit
1462!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
1463!! del vettore devono essere accordate. Il metodo controlla se il file è
1464!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1465SUBROUTINE geo_coord_read_unit(this, unit)
1466TYPE(geo_coord),INTENT(out) :: this !< oggetto da leggere
1467INTEGER, INTENT(in) :: unit !< unità da cui leggere
1468
1469CALL geo_coord_vect_read_unit((/this/), unit)
1470
1471END SUBROUTINE geo_coord_read_unit
1472
1473
1474!> Legge da un'unità di file il contenuto dell'oggetto \a this.
1475!! Il record da leggere deve essere stato scritto con la ::write_unit
1476!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
1477!! del vettore devono essere accordate. Il metodo controlla se il file è
1478!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1479SUBROUTINE geo_coord_vect_read_unit(this, unit)
1480TYPE(geo_coord) :: this(:) !< oggetto da leggere
1481INTEGER, INTENT(in) :: unit !< unità da cui leggere
1482
1483CHARACTER(len=40) :: form
1484INTEGER :: i
1485
1486INQUIRE(unit, form=form)
1487IF (form == 'FORMATTED') THEN
1488 read(unit,*) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1489!TODO bug gfortran compiler !
1490!missing values are unredeable when formatted
1491ELSE
1492 READ(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1493ENDIF
1494
1495END SUBROUTINE geo_coord_vect_read_unit
1496
1497
1498!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
1499!! Il record scritto potrà successivamente essere letto con la ::read_unit.
1500!! Il metodo controlla se il file è
1501!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1502SUBROUTINE geo_coord_write_unit(this, unit)
1503TYPE(geo_coord),INTENT(in) :: this !< oggetto da scrivere
1504INTEGER, INTENT(in) :: unit !< unità su cui scrivere
1505
1506CALL geo_coord_vect_write_unit((/this/), unit)
1507
1508END SUBROUTINE geo_coord_write_unit
1509
1510
1511!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
1512!! Il record scritto potrà successivamente essere letto con la ::read_unit.
1513!! Il metodo controlla se il file è
1514!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1515SUBROUTINE geo_coord_vect_write_unit(this, unit)
1516TYPE(geo_coord),INTENT(in) :: this(:) !< oggetto da scrivere
1517INTEGER, INTENT(in) :: unit !< unità su cui scrivere
1518
1519CHARACTER(len=40) :: form
1520INTEGER :: i
1521
1522INQUIRE(unit, form=form)
1523IF (form == 'FORMATTED') THEN
1524 WRITE(unit,*) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1525!TODO bug gfortran compiler !
1526!missing values are unredeable when formatted
1527ELSE
1528 WRITE(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
1529ENDIF
1530
1531END SUBROUTINE geo_coord_vect_write_unit
1532
1533
1534!> Restituisce la distanza in m tra 2 oggetti geo_coord.
1535!! La distanza è calcolata approssimativamente ed è valida per piccoli angoli.
1536FUNCTION geo_coord_dist(this, that) RESULT(dist)
1538TYPE(geo_coord), INTENT (IN) :: this !< primo punto
1539TYPE(geo_coord), INTENT (IN) :: that !< secondo punto
1540REAL(kind=fp_geo) :: dist !< distanza in metri
1541
1542REAL(kind=fp_geo) :: x,y
1543! Distanza approssimata, valida per piccoli angoli
1544
1545x = (getlon(this)-getlon(that))*cos(((getlat(this)+getlat(that))/2.)*degrad)
1546y = getlat(this)-getlat(that)
1547dist = sqrt(x**2 + y**2)*degrad*rearth
1548
1549END FUNCTION geo_coord_dist
1550
1551
1552!> Determina se il punto indicato da \a this è contenuto in un rettangolo.
1553!! Il rettangolo è orientato parallelamente agli assi del sistema,
1554!! i suoi vertici sud-ovest e nord-est sono specificati da altri due punti.
1555!! La funzione restituisce \c .TRUE. anche se il punto si trova sulla frontiera
1556!! del rettangolo.
1557!! Tutti gli oggetti devono essere già stati
1558!! convertiti ad un sistema di coordinate comune, altrimenti viene
1559!! restituito \c .FALSE. .
1560FUNCTION geo_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1561TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui determinare la posizione
1562TYPE(geo_coord),INTENT(IN) :: coordmin !< vertice sud-ovest del rettangolo
1563TYPE(geo_coord),INTENT(IN) :: coordmax !< vertice nord-est del rettangolo
1564LOGICAL :: res !< \c .TRUE. se \a this è dentro il rettangolo o sul bordo e \c .FALSE. se è fuori
1565
1566res = (geo_coord_ure(this, coordmin) .AND. geo_coord_lle(this,coordmax))
1567
1568END FUNCTION geo_coord_inside_rectang
1569
1570
1571! ===================
1572! == geo_coordvect ==
1573! ===================
1574!> Costruisce un oggetto \a geo_coordvect con i parametri opzionali forniti.
1575!! Se sono presenti \a lon e \a lat, inizializza le coordinate geografiche
1576!! ignorando \a utme e \a utmn, mentre se sono specificati \a utme e \a utmn
1577!! succede il contrario; non è possibile specificare le coordinate in entrambi
1578!! i sistemi, usare eventualmente \a to_geo. Se non viene passato nessun parametro
1579!! opzionale l'oggetto è inizializzato a valore mancante.
1580!! Il numero di punti dell'oggetto finale sarà uguale all'estensione
1581!! del più breve vettore della coppia fornita.
1582RECURSIVE SUBROUTINE geo_coordvect_init(this, lon, lat)
1583TYPE(geo_coordvect), INTENT(OUT) :: this !< oggetto da inizializzare
1584REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon(:) !< longitudine geografica
1585REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat(:) !< latitudine geografica
1586
1587IF (PRESENT(lon) .AND. PRESENT(lat)) THEN
1588 this%vsize = min(SIZE(lon), SIZE(lat))
1589 ALLOCATE(this%ll(this%vsize,2))
1590 this%ll(1:this%vsize,1) = lon(1:this%vsize)
1591 this%ll(1:this%vsize,2) = lat(1:this%vsize)
1592ELSE
1593 this%vsize = 0
1594 NULLIFY(this%ll)
1595ENDIF
1596this%vtype = 0 !?
1597
1598END SUBROUTINE geo_coordvect_init
1599
1600
1601!> Distrugge l'oggetto in maniera pulita, liberando l'eventuale spazio
1602!! dinamicamente allocato.
1603SUBROUTINE geo_coordvect_delete(this)
1604TYPE(geo_coordvect), INTENT(INOUT) :: this
1605
1606IF (ASSOCIATED(this%ll)) DEALLOCATE(this%ll)
1607this%vsize = 0
1608this%vtype = 0
1609
1610END SUBROUTINE geo_coordvect_delete
1611
1612
1613!> Restituisce il valore di uno o più componenti di un oggetto \a geo_coordvect.
1614!! Qualsiasi combinazione dei parametri opzionali è consentita; se
1615!! il tipo di coordinata richiesta non è stato inizializzato né calcolato,
1616!! restituisce il corrispondente valore mancante.
1617!! Se forniti, i parametri \a lon, \a lat, \a utme, \a utmn devono essere
1618!! dichiarati come puntatori che vengono
1619!! allocati dalla \a getval stessa e che devono poi essere
1620!! deallocati esplicitamente dal programma chiamante.
1621SUBROUTINE geo_coordvect_getval(this, lon, lat)
1622TYPE(geo_coordvect),INTENT(IN) :: this !< oggetto di cui restituire i componenti
1623REAL(kind=fp_geo), OPTIONAL, POINTER :: lon(:) !< longitudine geografica
1624REAL(kind=fp_geo), OPTIONAL, POINTER :: lat(:) !< latitudine geografica
1625
1626IF (PRESENT(lon)) THEN
1627 IF (ASSOCIATED(this%ll)) THEN
1628 ALLOCATE(lon(this%vsize))
1629 lon(:) = this%ll(1:this%vsize,1)
1630 ENDIF
1631ENDIF
1632IF (PRESENT(lat)) THEN
1633 IF (ASSOCIATED(this%ll)) THEN
1634 ALLOCATE(lat(this%vsize))
1635 lat(:) = this%ll(1:this%vsize,2)
1636 ENDIF
1637ENDIF
1638
1639END SUBROUTINE geo_coordvect_getval
1640
1641
1642SUBROUTINE geo_coordvect_import(this, unitsim, &
1643#ifdef HAVE_SHAPELIB
1644 shphandle, &
1645#endif
1646 nshp)
1647TYPE(geo_coordvect), INTENT(OUT) :: this
1648INTEGER,OPTIONAL,INTENT(IN) :: unitsim
1649#ifdef HAVE_SHAPELIB
1650TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
1651#endif
1652INTEGER,OPTIONAL,INTENT(IN) :: nshp
1653
1654REAL(kind=fp_geo),ALLOCATABLE :: llon(:), llat(:)
1655REAL(kind=fp_geo) :: lv1,lv2,lv3,lv4
1656INTEGER :: i, lvsize
1657CHARACTER(len=40) :: lname
1658#ifdef HAVE_SHAPELIB
1659TYPE(shpobject) :: shpobj
1660#endif
1661
1662IF (PRESENT(unitsim)) THEN
1663 ! Leggo l'intestazione
1664 READ(unitsim,*,END=10)lvsize,lv1,lv2,lv3,lv4,lname
1665 ALLOCATE(llon(lvsize+1), llat(lvsize+1))
1666 ! Leggo il poligono
1667 READ(unitsim,*)(llon(i),llat(i), i=1,lvsize)
1668 ! Lo chiudo se necessario
1669 IF (llon(1) /= llon(lvsize) .OR. llat(1) /= llat(lvsize)) THEN
1670 lvsize = lvsize + 1
1671 llon(lvsize) = llon(1)
1672 llat(lvsize) = llat(1)
1673 ENDIF
1674 ! Lo inserisco nel mio oggetto
1675 CALL init(this, lon=llon(1:lvsize), lat=llat(1:lvsize))
1676 this%vtype = geo_coordvect_polygon ! Sempre un poligono
1677
1678 DEALLOCATE(llon, llat)
1679 RETURN
168010 CALL raise_error('nella lettura del file '//trim(to_char(unitsim)))
1681 DEALLOCATE(llon, llat) ! End of file, ritorno un oggetto non assegnato
1682#ifdef HAVE_SHAPELIB
1683ELSE IF (PRESENT(shphandle) .AND. PRESENT(nshp)) THEN
1684 ! Leggo l'oggetto shape
1685 shpobj = shpreadobject(shphandle, nshp)
1686 IF (.NOT.shpisnull(shpobj)) THEN
1687 ! Lo inserisco nel mio oggetto
1688 CALL init(this, lon=real(shpobj%padfx,kind=fp_geo), &
1689 lat=real(shpobj%padfy,kind=fp_geo))
1690 this%vtype = shpobj%nshptype
1691 CALL shpdestroyobject(shpobj)
1692 ELSE
1693 CALL init(this)
1694 ENDIF
1695#endif
1696ENDIF
1697
1698END SUBROUTINE geo_coordvect_import
1699
1700
1701SUBROUTINE geo_coordvect_export(this, unitsim, &
1702#ifdef HAVE_SHAPELIB
1703 shphandle, &
1704#endif
1705 nshp)
1706TYPE(geo_coordvect), INTENT(INOUT) :: this
1707INTEGER,OPTIONAL,INTENT(IN) :: unitsim
1708#ifdef HAVE_SHAPELIB
1709TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
1710#endif
1711INTEGER,OPTIONAL,INTENT(IN) :: nshp
1712
1713INTEGER :: i, lnshp
1714#ifdef HAVE_SHAPELIB
1715TYPE(shpobject) :: shpobj
1716#endif
1717
1718IF (PRESENT(unitsim)) THEN
1719 IF (this%vsize > 0) THEN
1720 ! Scrivo l'intestazione
1721 WRITE(unitsim,*)SIZE(this%ll,1),-1.,5000.,-0.1,1.1,'Area'
1722 ! Scrivo il poligono
1723 WRITE(unitsim,*)(this%ll(i,1:2), i=1,this%vsize)
1724 ELSE
1725 CALL raise_warning('oggetto geo_coordvect vuoto, non scrivo niente in '// &
1726 trim(to_char(unitsim)))
1727 ENDIF
1728#ifdef HAVE_SHAPELIB
1729ELSE IF (PRESENT(shphandle)) THEN
1730 IF (PRESENT(nshp)) THEN
1731 lnshp = nshp
1732 ELSE
1733 lnshp = -1 ! -1 = append
1734 ENDIF
1735 ! Creo l'oggetto shape inizializzandolo con il mio oggetto
1736 shpobj = shpcreatesimpleobject(this%vtype, this%vsize, &
1737 REAL(this%ll(1:this%vsize,1),kind=fp_d), &
1738 REAL(this%ll(1:this%vsize,2),kind=fp_d))
1739 IF (.NOT.shpisnull(shpobj)) THEN
1740 ! Lo scrivo nello shapefile
1741 i=shpwriteobject(shphandle, lnshp, shpobj)
1742 CALL shpdestroyobject(shpobj)
1743 ENDIF
1744#endif
1745ENDIF
1746
1747END SUBROUTINE geo_coordvect_export
1748
1749!> Importa un vettore di oggetti \a geo_coordvect da un file in
1750!! formato testo o in formato \a shapefile.
1751!! Il parametro \a this è un puntatore che sarà allocato a cura del metodo stesso e
1752!! dovrà invece essere deallocato da parte del programma chiamante
1753!! dopo aver chiamato il metodo \a delete per ogni suo elemento.
1754!! In caso di errore nella fase iniziale di importazione,
1755!! \a this non verrà associato, e quindi è opportuno testare
1756!! \code
1757!! IF (ASSOCIATED(my_coord_vect)) THEN...
1758!! \endcode
1759!! nel programma chiamante
1760!! per intrappolare eventuale condizioni di errore (tipicamente file non
1761!! trovato o in un formato non compatibile).
1762!! Entrambi i formati di ingresso non contengono informazioni sul tipo
1763!! di coordinate dei dati
1764!! (per il formato shapefile è possibile solo con delle estensioni non standard),
1765!! per cui questa informazione, se desiderata, deve essere fornita dal
1766!! programma chiamante.
1767SUBROUTINE geo_coordvect_importvect(this, shpfilesim, shpfile)
1768TYPE(geo_coordvect),POINTER :: this(:) !< puntatore all'oggetto su cui importare i dati, viene allocato dalla \a import stessa
1769CHARACTER(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
1770CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile !< nome delllo shapefile, il parametro deve essere fornito solo se si vuole importare da un file di quel tipo
1771
1772REAL(kind=fp_geo) :: inu
1773REAL(kind=fp_d) :: minb(4), maxb(4)
1774INTEGER :: i, u, ns, lvsize, shptype, dbfnf, dbfnr
1775CHARACTER(len=40) :: lname
1776#ifdef HAVE_SHAPELIB
1777TYPE(shpfileobject) :: shphandle
1778#endif
1779
1780NULLIFY(this)
1781
1782IF (PRESENT(shpfilesim)) THEN
1783 u = getunit()
1784 OPEN(u, file=shpfilesim, status='old', err=30)
1785 ns = 0 ! Conto il numero di shape contenute
1786 DO WHILE(.true.)
1787 READ(u,*,END=10,ERR=20)lvsize,inu,inu,inu,inu,lname
1788 READ(u,*,END=20,ERR=20)(inu,inu, i=1,lvsize)
1789 ns = ns + 1
1790 ENDDO
179110 CONTINUE
1792 IF (ns > 0) THEN ! Alloco e leggo il mio oggetto
1793 ALLOCATE(this(ns))
1794 rewind(u)
1795 DO i = 1, ns
1796 CALL import(this(i), unitsim=u)
1797 ENDDO
1798 ENDIF
179920 CONTINUE
1800 CLOSE(u)
1801 IF (.NOT.ASSOCIATED(this)) THEN
1802 CALL raise_warning('file '//trim(shpfilesim)//' vuoto o corrotto')
1803 ENDIF
1804 RETURN
180530 CONTINUE
1806 CALL raise_error('Impossibile aprire il file '//trim(shpfile))
1807 RETURN
1808
1809ELSE IF (PRESENT(shpfile)) THEN
1810#ifdef HAVE_SHAPELIB
1811 shphandle = shpopen(trim(shpfile), 'rb')
1812 IF (shpfileisnull(shphandle)) THEN
1813 CALL raise_error('Impossibile aprire lo shapefile '//trim(shpfile))
1814 RETURN
1815 ENDIF
1816 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr) ! Ottengo le info sul file
1817 IF (ns > 0) THEN ! Alloco e leggo il mio oggetto
1818 ALLOCATE(this(ns))
1819 this(:)%vtype = shptype
1820 DO i = 1, ns
1821 CALL import(this(i), shphandle=shphandle, nshp=i-1)
1822 ENDDO
1823 ENDIF
1824 CALL shpclose(shphandle)
1825 RETURN
1826#endif
1827ENDIF
1828
1829END SUBROUTINE geo_coordvect_importvect
1830
1831
1832!> Esporta un vettore di oggetti \a geo_coordvect su un file in
1833!! formato testo o in formato \a shapefile.
1834SUBROUTINE geo_coordvect_exportvect(this, shpfilesim, shpfile, append)
1835TYPE(geo_coordvect) :: this(:) !< oggetto da esportare
1836CHARACTER(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
1837CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile !< nome dello shapefile, il parametro deve essere fornito solo se si vuole esportare su un file di quel tipo
1838!> sistema di coordinate (proiezione) dei dati,
1839!! usare i parametri \a ::proj_geo (default) o \a ::proj_utm,
1840!! ha senso se \a this, a seguito di una chiamata a \a ::to_geo o a \a ::to_utm,
1841!! contiene le coordinate in entrambi i sistemi, altrimenti i dati vengono
1842!! esportati automaticamente nel solo sistema disponibile
1843LOGICAL,INTENT(in),OPTIONAL :: append !< se è presente e vale \c .TRUE. , ::export accoda all'eventuale file esistente anziché sovrascriverlo
1844
1845REAL(kind=fp_d) :: minb(4), maxb(4)
1846INTEGER :: i, u, ns, shptype, dbfnf, dbfnr
1847LOGICAL :: lappend
1848#ifdef HAVE_SHAPELIB
1849TYPE(shpfileobject) :: shphandle
1850#endif
1851
1852IF (PRESENT(append)) THEN
1853 lappend = append
1854ELSE
1855 lappend = .false.
1856ENDIF
1857IF (PRESENT(shpfilesim)) THEN
1858 u = getunit()
1859 IF (lappend) THEN
1860 OPEN(u, file=shpfilesim, status='unknown', position='append', err=30)
1861 ELSE
1862 OPEN(u, file=shpfilesim, status='unknown', err=30)
1863 ENDIF
1864 DO i = 1, SIZE(this)
1865 CALL export(this(i), unitsim=u)
1866 ENDDO
1867 CLOSE(u)
1868 RETURN
186930 CONTINUE
1870 CALL raise_error('Impossibile aprire il file '//trim(shpfile))
1871 RETURN
1872ELSE IF (PRESENT(shpfile)) THEN
1873#ifdef HAVE_SHAPELIB
1874 IF (lappend) THEN
1875 shphandle = shpopen(trim(shpfile), 'r+b')
1876 IF (shpfileisnull(shphandle)) THEN ! shpopen funziona solo su file esistenti
1877 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
1878 ENDIF
1879 ELSE
1880 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
1881 ENDIF
1882 IF (shpfileisnull(shphandle)) THEN
1883 CALL raise_error('Impossibile aprire lo shapefile '//trim(shpfile))
1884 RETURN
1885 ENDIF
1886 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr) ! Ottengo le info sul file
1887 DO i = 1, SIZE(this)
1888 IF (i > ns .OR. lappend) THEN ! Append shape
1889 CALL export(this(i), shphandle=shphandle)
1890 ELSE ! Overwrite shape
1891 CALL export(this(i), shphandle=shphandle, nshp=i-1)
1892 ENDIF
1893 ENDDO
1894 CALL shpclose(shphandle)
1895 RETURN
1896#endif
1897ENDIF
1898
1899END SUBROUTINE geo_coordvect_exportvect
1900
1901
1902!> Determina se il punto indicato da \a this si trova
1903!! dentro o fuori dal poligono descritto dall'oggetto \a poly.
1904!! Funziona anche se la topologia di \a poly non è poligonale,
1905!! forzandone la chiusura; usa un algoritmo di ricerca del numero di
1906!! intersezioni, come indicato in
1907!! comp.graphics.algorithms FAQ (http://www.faqs.org/faqs/graphics/algorithms-faq/)
1908!! o in
1909!! http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
1910FUNCTION geo_coord_inside(this, poly) RESULT(inside)
1911TYPE(geo_coord), INTENT(IN) :: this !< oggetto di cui determinare la posizione
1912TYPE(geo_coordvect), INTENT(IN) :: poly !< poligono contenitore
1913LOGICAL :: inside
1914
1915INTEGER :: i, j, starti
1916
1917inside = .false.
1918IF (all(poly%ll(1,:) == poly%ll(poly%vsize,:))) THEN ! Poligono chiuso
1919 starti = 2
1920 j = 1
1921ELSE ! Poligono non chiuso
1922 starti = 1
1923 j = poly%vsize
1924ENDIF
1925DO i = starti, poly%vsize
1926 IF ((poly%ll(i,2) <= getlat(this) .AND. &
1927 getlat(this) < poly%ll(j,2)) .OR. &
1928 (poly%ll(j,2) <= getlat(this) .AND. &
1929 getlat(this) < poly%ll(i,2))) THEN
1930 IF (getlon(this) < (poly%ll(j,1) - poly%ll(i,1)) * &
1931 (getlat(this) - poly%ll(i,2)) / &
1932 (poly%ll(j,2) - poly%ll(i,2)) + poly%ll(i,1)) THEN
1933 inside = .NOT. inside
1934 ENDIF
1935 ENDIF
1936 j = i
1937ENDDO
1938
1939END FUNCTION geo_coord_inside
1940
1941
1942ELEMENTAL FUNCTION c_e_geo_coord(this) result (res)
1943TYPE(geo_coord),INTENT(in) :: this
1944LOGICAL :: res
1945
1946res = .not. this == geo_coord_miss
1947
1948end FUNCTION c_e_geo_coord
1949
1950
1951character(len=80) function to_char_geo_coord(this)
1952TYPE(geo_coord),INTENT(in) :: this
1953
1954to_char_geo_coord = "GEO_COORD: Lon="// &
1955 trim(to_char(getlon(this),miss="Missing lon",form="(f11.5)"))//&
1956 " Lat="// &
1957 trim(to_char(getlat(this),miss="Missing lat",form="(f11.5)"))
1958
1959end function to_char_geo_coord
1960
1961
1962subroutine display_geo_coord(this)
1963TYPE(geo_coord),INTENT(in) :: this
1964
1965print*,trim(to_char(this))
1966
1967end subroutine display_geo_coord
1968
1969
1970END MODULE geo_coord_class
Detructors for the two classes.
Export one or more geo_coordvect objects to a plain text file or to a file in ESRI/Shapefile format.
Methods for returning the value of object members.
Import one or more geo_coordvect objects from a plain text file or for a file in ESRI/Shapefile forma...
Constructors for the two classes.
Determine whether a point lies inside a polygon or a rectangle.
Read a single geo_coord object or an array of geo_coord objects from a Fortran FORMATTED or UNFORMATT...
Represent geo_coord object in a pretty string.
Write a single geo_coord object or an array of geo_coord objects to a Fortran FORMATTED or UNFORMATTE...
Utilities for CHARACTER variables.
Costanti fisiche (DOUBLEPRECISION).
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.