libsim Versione 7.2.6

◆ georef_coord_read_unit()

subroutine georef_coord_read_unit ( type(georef_coord), intent(out) this,
integer, intent(in) unit )

Legge da un'unità di file il contenuto dell'oggetto this.

Il record da leggere deve essere stato scritto con la write_unit e, nel caso this sia un vettore, la lunghezza del record e quella del vettore devono essere accordate. Il metodo controlla se il file è aperto per un I/O formattato o non formattato e fa la cosa giusta.

Parametri
[out]thisoggetto da leggere
[in]unitunità da cui leggere

Definizione alla linea 850 del file georef_coord_class.F90.

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