libsim Versione 7.2.6

◆ georef_coord_vect_write_unit()

subroutine georef_coord_vect_write_unit ( type(georef_coord), dimension(:), 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 900 del file georef_coord_class.F90.

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