libsim Versione 7.2.6

◆ georef_coord_write_unit()

subroutine georef_coord_write_unit ( type(georef_coord), intent(in) this,
integer, intent(in) unit )
private

Scrive su un'unità di file il contenuto dell'oggetto this.

Il record scritto potrà successivamente essere letto con la read_unit. Il metodo controlla se il file è aperto per un I/O formattato o non formattato e fa la cosa giusta.

Parametri
[in]thisoggetto da scrivere
[in]unitunità su cui scrivere

Definizione alla linea 887 del file georef_coord_class.F90.

888! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
889! authors:
890! Davide Cesari <dcesari@arpa.emr.it>
891! Paolo Patruno <ppatruno@arpa.emr.it>
892
893! This program is free software; you can redistribute it and/or
894! modify it under the terms of the GNU General Public License as
895! published by the Free Software Foundation; either version 2 of
896! the License, or (at your option) any later version.
897
898! This program is distributed in the hope that it will be useful,
899! but WITHOUT ANY WARRANTY; without even the implied warranty of
900! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
901! GNU General Public License for more details.
902
903! You should have received a copy of the GNU General Public License
904! along with this program. If not, see <http://www.gnu.org/licenses/>.
905#include "config.h"
906!> This module defines objects describing georeferenced sparse points
907!! possibly with topology and projection information. This module
908!! defines two classes, \a georef_coord, which represents a single
909!! georeferenced point on the Earth, and \a georef_coord_array which
910!! defines a set of points with a topological relation.
911!!
912!! Both classes have \a PRIVATE members, so that they cannot be
913!! manipulated directly, but only through the proper methods.
914!!
915!! It is also possible to dafine a dynamically extendible array of \a
916!! georef_coord_array objects, of type \a arrayof_georef_coord_array,
917!! suitable for importing/exporting data from/to a shapefile.
918!!
919!! \ingroup base
921USE err_handling
924USE geo_proj_class
925#ifdef HAVE_SHAPELIB
926USE shapelib
927#endif
928IMPLICIT NONE
929
930!> Derive type defining a single georeferenced point, either in
931!! geodetic or in projected coordinates. The object has no information
932!! about the georeferenced coordinate system associated, it must be
933!! kept separately by the user.
934TYPE georef_coord
935 PRIVATE
936 DOUBLE PRECISION :: x=dmiss, y=dmiss
937END TYPE georef_coord
938
939!> Missing value for georef_coord
940TYPE(georef_coord),PARAMETER :: georef_coord_miss=georef_coord(dmiss,dmiss)
941
942!> Derived type defining a one-dimensional array of georeferenced points
943!! with an associated topology (isolated point, arc, polygon, group of
944!! points), possibly broken into parts and with an associated
945!! georeferenced coordinate system.
947 PRIVATE
948 INTEGER,ALLOCATABLE :: parts(:)
949 TYPE(georef_coord),ALLOCATABLE :: coord(:)
950 INTEGER :: topo=imiss
951 TYPE(geo_proj) :: proj
952 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
953 LOGICAL :: bbox_updated=.false.
954END TYPE georef_coord_array
955
956INTEGER,PARAMETER :: georef_coord_array_point = 1 !< Topology for georef_coord_array (from shapelib): isolated point
957INTEGER,PARAMETER :: georef_coord_array_arc = 3 !< Topology for georef_coord_array (from shapelib): arc (multiple arcs unsupported)
958INTEGER,PARAMETER :: georef_coord_array_polygon = 5 !< Topology for georef_coord_array (from shapelib): polygon (necessarily closed, multiple polygons unsupported)
959INTEGER,PARAMETER :: georef_coord_array_multipoint = 8 !< Topology for georef_coord_array (from shapelib): group of points
960
961
962!> Detructors for the two classes.
963!! They clean up all the information associated with the corresponding
964!! objects.
965INTERFACE delete
966 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
967END INTERFACE
968
969!> Check missing value.
970INTERFACE c_e
971 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
972END INTERFACE
973
974!> Methods for returning the value of object members.
975INTERFACE getval
976 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
977END INTERFACE
978
979INTERFACE compute_bbox
980 MODULE PROCEDURE georef_coord_array_compute_bbox
981END INTERFACE
982
983!> Logical equality operator.
984INTERFACE OPERATOR (==)
985 MODULE PROCEDURE georef_coord_eq
986END INTERFACE
987
988!> Logical inequality operator.
989INTERFACE OPERATOR (/=)
990 MODULE PROCEDURE georef_coord_ne
991END INTERFACE
992
993!> Logical greater-equal operator. Returns true if the first point
994!! lies to the northeast of the second
995INTERFACE OPERATOR (>=)
996 MODULE PROCEDURE georef_coord_ge
997END INTERFACE
998
999!> Logical less-equal operator. Returns true if the first point
1000!! lies to the southwest of the second
1001INTERFACE OPERATOR (<=)
1002 MODULE PROCEDURE georef_coord_le
1003END INTERFACE
1004
1005#ifdef HAVE_SHAPELIB
1006!> Import an array of \a georef_coord_array objects from a file
1007!! in ESRI/Shapefile format.
1008INTERFACE import
1009 MODULE PROCEDURE arrayof_georef_coord_array_import
1010END INTERFACE
1011
1012!> Export an array of \a georef_coord_array objects to a file
1013!! in ESRI/Shapefile format.
1014INTERFACE export
1015 MODULE PROCEDURE arrayof_georef_coord_array_export
1016END INTERFACE
1017#endif
1018
1019!> Read a single \a georef_coord object or an array of \a georef_coord objects
1020!! from a Fortran \c FORMATTED or \c UNFORMATTED file.
1021INTERFACE read_unit
1022 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1023END INTERFACE
1024
1025!> Write a single \a georef_coord object or an array of \a georef_coord objects
1026!! to a Fortran \c FORMATTED or \c UNFORMATTED file.
1027INTERFACE write_unit
1028 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1029END INTERFACE
1030
1031!> Determine whether a point lies inside a polygon or a rectangle.
1032INTERFACE inside
1033 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1034END INTERFACE
1035
1036!> Compute the distance in m between two points
1037INTERFACE dist
1038 MODULE PROCEDURE georef_coord_dist
1039END INTERFACE
1040
1041#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1042#define ARRAYOF_TYPE arrayof_georef_coord_array
1043!define ARRAYOF_ORIGEQ 0
1044#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1045#include "arrayof_pre.F90"
1046! from arrayof
1047PUBLIC insert, append, remove, packarray
1048
1049PRIVATE
1050PUBLIC georef_coord, georef_coord_miss, &
1051 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1052 georef_coord_array_polygon, georef_coord_array_multipoint, &
1053 delete, c_e, getval, compute_bbox, OPERATOR(==), OPERATOR(/=), OPERATOR(>=), OPERATOR(<=), &
1054#ifdef HAVE_SHAPELIB
1055 import, export, &
1056#endif
1058 georef_coord_new, georef_coord_array_new
1059
1060CONTAINS
1061
1062#include "arrayof_post.F90"
1063
1064! ===================
1065! == georef_coord ==
1066! ===================
1067!> Construct a \a georef_coord object with the optional parameters provided.
1068!! If coordinates are not provided the object obtained is empty
1069!! (missing, see c_e function).
1070FUNCTION georef_coord_new(x, y) RESULT(this)
1071DOUBLE PRECISION,INTENT(in),OPTIONAL :: x !< x coordinate
1072DOUBLE PRECISION,INTENT(in),OPTIONAL :: y !< y coordinate
1073TYPE(georef_coord) :: this
1074
1075CALL optio(x, this%x)
1076CALL optio(y, this%y)
1077
1078END FUNCTION georef_coord_new
1079
1080
1081SUBROUTINE georef_coord_delete(this)
1082TYPE(georef_coord),INTENT(inout) :: this
1083
1084this%x = dmiss
1085this%y = dmiss
1086
1087END SUBROUTINE georef_coord_delete
1088
1089
1090ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1091TYPE(georef_coord),INTENT(in) :: this
1092LOGICAL :: res
1093
1094res = .NOT. this == georef_coord_miss
1095
1096END FUNCTION georef_coord_c_e
1097
1098
1099!> Query a \a georef_coord object.
1100!! This is the correct way to retrieve the contents of a \a
1101!! georef_coord object, since its members are declared as \a
1102!! PRIVATE. It is declared as \a ELEMENTAL, thus it works also on
1103!! arrays of any shape and, in that case, the result will hae the same
1104!! shape as \a this.
1105ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1106TYPE(georef_coord),INTENT(in) :: this !< object to query
1107DOUBLE PRECISION,INTENT(out),OPTIONAL :: x !< x-coordinate
1108DOUBLE PRECISION,INTENT(out),OPTIONAL :: y !< y-coordinate
1109
1110IF (PRESENT(x)) x = this%x
1111IF (PRESENT(y)) y = this%y
1112
1113END SUBROUTINE georef_coord_getval
1114
1115
1116!> Query a \a georef_coord object associating a geographical projection to it.
1117!! This method allow to interpret the \a x,y coordinates of a \a
1118!! georef_coord object as projected on a specified geographical
1119!! projection system and retrieve the geodetic longitude and/or
1120!! latitude associated to them. When \a x or \y are requested it works
1121!! as the basic get_val method. It is declared as \a ELEMENTAL, thus
1122!! it works also on arrays of any shape and, in that case, the result
1123!! will hae the same shape as \a this.
1124ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1125TYPE(georef_coord),INTENT(in) :: this !< object to query
1126TYPE(geo_proj),INTENT(in) :: proj !< geographical projection associated to coordinates of \a this
1127DOUBLE PRECISION,INTENT(out),OPTIONAL :: x !< x-coordinate
1128DOUBLE PRECISION,INTENT(out),OPTIONAL :: y !< y-coordinate
1129DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon !< geodetic longitude
1130DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat !< geodetic latitude
1131
1132DOUBLE PRECISION :: llon, llat
1133
1134IF (PRESENT(x)) x = this%x
1135IF (PRESENT(y)) y = this%y
1136IF (PRESENT(lon) .OR. present(lat)) THEN
1137 CALL unproj(proj, this%x, this%y, llon, llat)
1138 IF (PRESENT(lon)) lon = llon
1139 IF (PRESENT(lat)) lat = llat
1140ENDIF
1141
1142END SUBROUTINE georef_coord_proj_getval
1143
1144
1145! document and improve
1146ELEMENTAL FUNCTION getlat(this)
1147TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1148DOUBLE PRECISION :: getlat ! latitudine geografica
1149
1150getlat = this%y ! change!!!
1151
1152END FUNCTION getlat
1153
1154! document and improve
1155ELEMENTAL FUNCTION getlon(this)
1156TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1157DOUBLE PRECISION :: getlon ! longitudine geografica
1158
1159getlon = this%x ! change!!!
1160
1161END FUNCTION getlon
1162
1163
1164ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1165TYPE(georef_coord),INTENT(IN) :: this, that
1166LOGICAL :: res
1167
1168res = (this%x == that%x .AND. this%y == that%y)
1169
1170END FUNCTION georef_coord_eq
1171
1172
1173ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1174TYPE(georef_coord),INTENT(IN) :: this, that
1175LOGICAL :: res
1176
1177res = (this%x >= that%x .AND. this%y >= that%y)
1178
1179END FUNCTION georef_coord_ge
1180
1181
1182ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1183TYPE(georef_coord),INTENT(IN) :: this, that
1184LOGICAL :: res
1185
1186res = (this%x <= that%x .AND. this%y <= that%y)
1187
1188END FUNCTION georef_coord_le
1189
1190
1191ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1192TYPE(georef_coord),INTENT(IN) :: this, that
1193LOGICAL :: res
1194
1195res = .NOT.(this == that)
1196
1197END FUNCTION georef_coord_ne
1198
1199
1200!> Legge da un'unità di file il contenuto dell'oggetto \a this.
1201!! Il record da leggere deve essere stato scritto con la ::write_unit
1202!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
1203!! del vettore devono essere accordate. Il metodo controlla se il file è
1204!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1205SUBROUTINE georef_coord_read_unit(this, unit)
1206TYPE(georef_coord),INTENT(out) :: this !< oggetto da leggere
1207INTEGER, INTENT(in) :: unit !< unità da cui leggere
1208
1209CALL georef_coord_vect_read_unit((/this/), unit)
1210
1211END SUBROUTINE georef_coord_read_unit
1212
1213
1214!> Legge da un'unità di file il contenuto dell'oggetto \a this.
1215!! Il record da leggere deve essere stato scritto con la ::write_unit
1216!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
1217!! del vettore devono essere accordate. Il metodo controlla se il file è
1218!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1219SUBROUTINE georef_coord_vect_read_unit(this, unit)
1220TYPE(georef_coord) :: this(:) !< oggetto da leggere
1221INTEGER, INTENT(in) :: unit !< unità da cui leggere
1222
1223CHARACTER(len=40) :: form
1224INTEGER :: i
1225
1226INQUIRE(unit, form=form)
1227IF (form == 'FORMATTED') THEN
1228 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1229!TODO bug gfortran compiler !
1230!missing values are unredeable when formatted
1231ELSE
1232 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1233ENDIF
1234
1235END SUBROUTINE georef_coord_vect_read_unit
1236
1237
1238!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
1239!! Il record scritto potrà successivamente essere letto con la ::read_unit.
1240!! Il metodo controlla se il file è
1241!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1242SUBROUTINE georef_coord_write_unit(this, unit)
1243TYPE(georef_coord),INTENT(in) :: this !< oggetto da scrivere
1244INTEGER, INTENT(in) :: unit !< unità su cui scrivere
1245
1246CALL georef_coord_vect_write_unit((/this/), unit)
1247
1248END SUBROUTINE georef_coord_write_unit
1249
1250
1251!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
1252!! Il record scritto potrà successivamente essere letto con la ::read_unit.
1253!! Il metodo controlla se il file è
1254!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
1255SUBROUTINE georef_coord_vect_write_unit(this, unit)
1256TYPE(georef_coord),INTENT(in) :: this(:) !< oggetto da scrivere
1257INTEGER, INTENT(in) :: unit !< unità su cui scrivere
1258
1259CHARACTER(len=40) :: form
1260INTEGER :: i
1261
1262INQUIRE(unit, form=form)
1263IF (form == 'FORMATTED') THEN
1264 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1265!TODO bug gfortran compiler !
1266!missing values are unredeable when formatted
1267ELSE
1268 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1269ENDIF
1270
1271END SUBROUTINE georef_coord_vect_write_unit
1272
1273
1274!> Restituisce la distanza in m tra 2 oggetti georef_coord.
1275!! La distanza è calcolata approssimativamente ed è valida per piccoli angoli.
1276FUNCTION georef_coord_dist(this, that) RESULT(dist)
1278TYPE(georef_coord), INTENT (IN) :: this !< primo punto
1279TYPE(georef_coord), INTENT (IN) :: that !< secondo punto
1280DOUBLE PRECISION :: dist !< distanza in metri
1281
1282DOUBLE PRECISION :: x,y
1283! Distanza approssimata, valida per piccoli angoli
1284
1285x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1286y = (this%y-that%y)
1287dist = sqrt(x**2 + y**2)*degrad*rearth
1288
1289END FUNCTION georef_coord_dist
1290
1291
1292!> Determines whether the point \a this lies inside a specified rectangle.
1293!! The rectangle is oriented parallely to the coordinate system, its
1294!! lower-left and upper-right vertices are specified by the two
1295!! arguments. The function returns \c .TRUE. also in the case it lies
1296!! exactly on the border of the rectangle.
1297FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1298TYPE(georef_coord),INTENT(IN) :: this !< point to test
1299TYPE(georef_coord),INTENT(IN) :: coordmin !< lower-left vertex
1300TYPE(georef_coord),INTENT(IN) :: coordmax !< upper-right vertex
1301LOGICAL :: res
1302
1303res = (this >= coordmin .AND. this <= coordmax)
1304
1305END FUNCTION georef_coord_inside_rectang
1306
1307
1308! ========================
1309! == georef_coord_array ==
1310! ========================
1311!> Construct a \a georef_coord_array object with the optional parameters provided.
1312!! If coordinates are not provided the object obtained is empty
1313!! (missing, see c_e function), if coordinate arrays are of different
1314!! lengths the \a georef_coord_array is initialised to the shortest
1315!! length provided.
1316FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1317DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:) !< x coordinate array
1318DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:) !< y coordinate array
1319INTEGER,INTENT(in),OPTIONAL :: topo !< topology of the object, one of the \a georef_coord_array_* constants
1320TYPE(geo_proj),INTENT(in),OPTIONAL :: proj !< geographical projection associated to the coordinates
1321TYPE(georef_coord_array) :: this
1322
1323INTEGER :: lsize
1324
1325IF (PRESENT(x) .AND. PRESENT(y)) THEN
1326 lsize = min(SIZE(x), SIZE(y))
1327 ALLOCATE(this%coord(lsize))
1328 this%coord(1:lsize)%x = x(1:lsize)
1329 this%coord(1:lsize)%y = y(1:lsize)
1330ENDIF
1331this%topo = optio_l(topo)
1332IF (PRESENT(proj)) this%proj = proj
1333
1334END FUNCTION georef_coord_array_new
1335
1336
1337SUBROUTINE georef_coord_array_delete(this)
1338TYPE(georef_coord_array),INTENT(inout) :: this
1339
1340TYPE(georef_coord_array) :: lobj
1341
1342this = lobj
1343
1344END SUBROUTINE georef_coord_array_delete
1345
1346
1347ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1348TYPE(georef_coord_array),INTENT(in) :: this
1349LOGICAL :: res
1350
1351res = ALLOCATED(this%coord)
1352
1353END FUNCTION georef_coord_array_c_e
1354
1355
1356!> Query a \a georef_coord_array object.
1357!! This is the correct way to retrieve the contents of a \a
1358!! georef_coord_array object, since its members are declared as \a
1359!! PRIVATE.
1360SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1361TYPE(georef_coord_array),INTENT(in) :: this !< object to query
1362DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:) !< x-coordinate
1363DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:) !< y-coordinate
1364! allocatable per vedere di nascosto l'effetto che fa
1365INTEGER,OPTIONAL,INTENT(out) :: topo !< topology associated with the coordinates
1366TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj !< geographical projection
1367
1368
1369IF (PRESENT(x)) THEN
1370 IF (ALLOCATED(this%coord)) THEN
1371 x = this%coord%x
1372 ENDIF
1373ENDIF
1374IF (PRESENT(y)) THEN
1375 IF (ALLOCATED(this%coord)) THEN
1376 y = this%coord%y
1377 ENDIF
1378ENDIF
1379IF (PRESENT(topo)) topo = this%topo
1380IF (PRESENT(proj)) proj = this%proj ! warning proj has no missing value yet
1381
1382END SUBROUTINE georef_coord_array_getval
1383
1384
1385!> Compute the bounding box of each shape in \a georef_coord_array object.
1386!! The bounding box is computed and stored in the object, it is used
1387!! by the inside() function for speedup; after it is computed the
1388!! object cannot be changed, otherwise the bounding box will not be
1389!! valid.
1390SUBROUTINE georef_coord_array_compute_bbox(this)
1391TYPE(georef_coord_array),INTENT(inout) :: this !< object to manipulate
1392
1393IF (ALLOCATED(this%coord)) THEN
1394 this%bbox(1)%x = minval(this%coord(:)%x)
1395 this%bbox(1)%y = minval(this%coord(:)%y)
1396 this%bbox(2)%x = maxval(this%coord(:)%x)
1397 this%bbox(2)%y = maxval(this%coord(:)%y)
1398 this%bbox_updated = .true.
1399ENDIF
1400
1401END SUBROUTINE georef_coord_array_compute_bbox
1402
1403#ifdef HAVE_SHAPELIB
1404! internal method for importing a single shape
1405SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1406TYPE(georef_coord_array),INTENT(OUT) :: this
1407TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1408INTEGER,INTENT(IN) :: nshp
1409
1410TYPE(shpobject) :: shpobj
1411
1412! read shape object
1413shpobj = shpreadobject(shphandle, nshp)
1414IF (.NOT.shpisnull(shpobj)) THEN
1415! import it in georef_coord object
1416 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1417 topo=shpobj%nshptype)
1418 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1419 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1420 ELSE IF (ALLOCATED(this%parts)) THEN
1421 DEALLOCATE(this%parts)
1422 ENDIF
1423 CALL shpdestroyobject(shpobj)
1424 CALL compute_bbox(this)
1425ENDIF
1426
1427
1428END SUBROUTINE georef_coord_array_import
1429
1430
1431! internal method for exporting a single shape
1432SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1433TYPE(georef_coord_array),INTENT(in) :: this
1434TYPE(shpfileobject),INTENT(inout) :: shphandle
1435INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1436
1437INTEGER :: i
1438TYPE(shpobject) :: shpobj
1439
1440IF (ALLOCATED(this%coord)) THEN
1441 IF (ALLOCATED(this%parts)) THEN
1442 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1443 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1444 ELSE
1445 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1446 this%coord(:)%x, this%coord(:)%y)
1447 ENDIF
1448ELSE
1449 RETURN
1450ENDIF
1451
1452IF (.NOT.shpisnull(shpobj)) THEN
1453 i = shpwriteobject(shphandle, nshp, shpobj)
1454 CALL shpdestroyobject(shpobj)
1455ENDIF
1456
1457END SUBROUTINE georef_coord_array_export
1458
1459
1460!> Import an array of \a georef_coord_array objects from a file
1461!! in ESRI/Shapefile format. The \a this argument is an uninitialised
1462!! \a arrayof_georef_coord_array, every element of which,
1463!! this%array(n), is of type \a georef_coord_array and, on return,
1464!! will contain information from the n-th shape of the file. Topology
1465!! information and possible polygon parts are imported as well, while
1466!! no projection information, even if available, is imported. An
1467!! error condition while opening the file can be detected by checking
1468!! .NOT.ASSOCIATED(this%array), while an error reading shape \a n can
1469!! be detected by checking .NOT.c_e(this%array(n)).
1470SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1471TYPE(arrayof_georef_coord_array),INTENT(out) :: this !< uninitialised array object
1472CHARACTER(len=*),INTENT(in) :: shpfile !< name of shapefile (with or without extension)
1473
1474REAL(kind=fp_d) :: minb(4), maxb(4)
1475INTEGER :: i, ns, shptype, dbfnf, dbfnr
1476TYPE(shpfileobject) :: shphandle
1477
1478shphandle = shpopen(trim(shpfile), 'rb')
1479IF (shpfileisnull(shphandle)) THEN
1480 ! log here
1481 CALL raise_error()
1482 RETURN
1483ENDIF
1484
1485! get info about file
1486CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1487IF (ns > 0) THEN ! allocate and read the object
1488 CALL insert(this, nelem=ns)
1489 DO i = 1, ns
1490 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1491 ENDDO
1492ENDIF
1493
1494CALL shpclose(shphandle)
1495! pack object to save memory
1496CALL packarray(this)
1497
1498END SUBROUTINE arrayof_georef_coord_array_import
1499
1500
1501!> Export an array of \a georef_coord_array objects to a file
1502!! in ESRI/Shapefile format. All the \a this%arraysize shapes
1503!! contained in \a this are exported to the requested
1504!! shapefile. Topology information and possible polygon parts are
1505!! exported as well, while projection information is ignored.
1506SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1507TYPE(arrayof_georef_coord_array),INTENT(in) :: this !< array object to be exported
1508CHARACTER(len=*),INTENT(in) :: shpfile !< name of shapefile (with or without extension)
1509
1510INTEGER :: i
1511TYPE(shpfileobject) :: shphandle
1512
1513IF (this%arraysize > 0) THEN
1514 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1515ELSE
1516 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1517ENDIF
1518IF (shpfileisnull(shphandle)) THEN
1519 ! log here
1520 CALL raise_error()
1521 RETURN
1522ENDIF
1523
1524DO i = 1, this%arraysize
1525 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1526ENDDO
1527
1528CALL shpclose(shphandle)
1529
1530END SUBROUTINE arrayof_georef_coord_array_export
1531#endif
1532
1533!> Determines whether the point \a this lies inside the polygon \a poly.
1534!! The polygon is forced to be closed if it is not already the case,
1535!! and there is no check about the topology of \a poly to really be of
1536!! polygon type. It works also with polygons in parts (as from
1537!! shapefile specification) defining either multiple polygons or
1538!! polygons with holes.
1539!!
1540!! The method used consists in counting the number of intersections as
1541!! indicated in comp.graphics.algorithms FAQ
1542!! (http://www.faqs.org/faqs/graphics/algorithms-faq/) or in
1543!! http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html.
1544FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1545TYPE(georef_coord), INTENT(IN) :: this !< oggetto di cui determinare la posizione
1546TYPE(georef_coord_array), INTENT(IN) :: poly !< poligono contenitore
1547LOGICAL :: inside
1548
1549INTEGER :: i
1550
1551inside = .false.
1552IF (.NOT.c_e(this)) RETURN
1553IF (.NOT.ALLOCATED(poly%coord)) RETURN
1554! if outside bounding box stop here
1555IF (poly%bbox_updated) THEN
1556 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1557ENDIF
1558
1559IF (ALLOCATED(poly%parts)) THEN
1560 DO i = 1, SIZE(poly%parts)-1
1561 inside = inside .NEQV. pointinpoly(this%x, this%y, &
1562 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1563 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1564 ENDDO
1565 IF (SIZE(poly%parts) > 0) THEN ! safety check
1566 inside = inside .NEQV. pointinpoly(this%x, this%y, &
1567 poly%coord(poly%parts(i)+1:)%x, &
1568 poly%coord(poly%parts(i)+1:)%y)
1569 ENDIF
1570
1571ELSE
1572 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1573 inside = pointinpoly(this%x, this%y, &
1574 poly%coord(:)%x, poly%coord(:)%y)
1575ENDIF
1576
1577CONTAINS
1578
1579FUNCTION pointinpoly(x, y, px, py)
1580DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1581LOGICAL :: pointinpoly
1582
1583INTEGER :: i, j, starti
1584
1585pointinpoly = .false.
1586
1587IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1588 starti = 2
1589 j = 1
1590ELSE ! unclosed polygon
1591 starti = 1
1592 j = SIZE(px)
1593ENDIF
1594DO i = starti, SIZE(px)
1595 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1596 (py(j) <= y .AND. y < py(i))) THEN
1597 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1598 pointinpoly = .NOT. pointinpoly
1599 ENDIF
1600 ENDIF
1601 j = i
1602ENDDO
1603
1604END FUNCTION pointinpoly
1605
1606END FUNCTION georef_coord_inside
1607
1608
1609
1610END MODULE georef_coord_class
Compute forward coordinate transformation from geographical system to projected system.
Compute backward coordinate transformation from projected system to geographical system.
Quick method to append an element to the array.
Detructors for the two classes.
Compute the distance in m between two points.
Export an array of georef_coord_array objects to a file in ESRI/Shapefile format.
Methods for returning the value of object members.
Import an array of georef_coord_array objects from a file in ESRI/Shapefile format.
Method for inserting elements of the array at a desired position.
Determine whether a point lies inside a polygon or a rectangle.
Method for packing the array object reducing at a minimum the memory occupation, without destroying i...
Read a single georef_coord object or an array of georef_coord objects from a Fortran FORMATTED or UNF...
Method for removing elements of the array at a desired position.
Write a single georef_coord object or an array of georef_coord objects to a Fortran FORMATTED or UNFO...
Generic subroutine for checking OPTIONAL parameters.
Costanti fisiche (DOUBLEPRECISION).
Gestione degli errori.
This module defines objects describing georeferenced sparse points possibly with topology and project...
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Derived type defining a one-dimensional array of georeferenced points with an associated topology (is...
Derive type defining a single georeferenced point, either in geodetic or in projected coordinates.

Generated with Doxygen.