libsim Versione 7.2.6

◆ georef_coord_proj_getval()

elemental subroutine georef_coord_proj_getval ( type(georef_coord), intent(in) this,
type(geo_proj), intent(in) proj,
double precision, intent(out), optional x,
double precision, intent(out), optional y,
double precision, intent(out), optional lon,
double precision, intent(out), optional lat )

Query a georef_coord object associating a geographical projection to it.

This method allow to interpret the x,y coordinates of a georef_coord object as projected on a specified geographical projection system and retrieve the geodetic longitude and/or latitude associated to them. When x or \y are requested it works as the basic get_val method. 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.

Parametri
[in]thisobject to query
[in]projgeographical projection associated to coordinates of this
[out]xx-coordinate
[out]yy-coordinate
[out]longeodetic longitude
[out]latgeodetic latitude

Definizione alla linea 769 del file georef_coord_class.F90.

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