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