libsim Versione 7.1.11
|
◆ georef_coord_proj_getval()
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.
Definizione alla linea 775 del file georef_coord_class.F90. 776! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
777! authors:
778! Davide Cesari <dcesari@arpa.emr.it>
779! Paolo Patruno <ppatruno@arpa.emr.it>
780
781! This program is free software; you can redistribute it and/or
782! modify it under the terms of the GNU General Public License as
783! published by the Free Software Foundation; either version 2 of
784! the License, or (at your option) any later version.
785
786! This program is distributed in the hope that it will be useful,
787! but WITHOUT ANY WARRANTY; without even the implied warranty of
788! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
789! GNU General Public License for more details.
790
791! You should have received a copy of the GNU General Public License
792! along with this program. If not, see <http://www.gnu.org/licenses/>.
793#include "config.h"
794
812USE geo_proj_class
813#ifdef HAVE_SHAPELIB
814USE shapelib
815#endif
816IMPLICIT NONE
817
823 PRIVATE
824 DOUBLE PRECISION :: x=dmiss, y=dmiss
826
829
835 PRIVATE
836 INTEGER,ALLOCATABLE :: parts(:)
837 TYPE(georef_coord),ALLOCATABLE :: coord(:)
838 INTEGER :: topo=imiss
839 TYPE(geo_proj) :: proj
840 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
841 LOGICAL :: bbox_updated=.false.
843
844INTEGER,PARAMETER :: georef_coord_array_point = 1
845INTEGER,PARAMETER :: georef_coord_array_arc = 3
846INTEGER,PARAMETER :: georef_coord_array_polygon = 5
847INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
848
849
854 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
855END INTERFACE
856
859 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
860END INTERFACE
861
864 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
865END INTERFACE
866
867INTERFACE compute_bbox
868 MODULE PROCEDURE georef_coord_array_compute_bbox
869END INTERFACE
870
872INTERFACE OPERATOR (==)
873 MODULE PROCEDURE georef_coord_eq
874END INTERFACE
875
877INTERFACE OPERATOR (/=)
878 MODULE PROCEDURE georef_coord_ne
879END INTERFACE
880
883INTERFACE OPERATOR (>=)
884 MODULE PROCEDURE georef_coord_ge
885END INTERFACE
886
889INTERFACE OPERATOR (<=)
890 MODULE PROCEDURE georef_coord_le
891END INTERFACE
892
893#ifdef HAVE_SHAPELIB
894
896INTERFACE import
897 MODULE PROCEDURE arrayof_georef_coord_array_import
898END INTERFACE
899
903 MODULE PROCEDURE arrayof_georef_coord_array_export
904END INTERFACE
905#endif
906
910 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
911END INTERFACE
912
916 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
917END INTERFACE
918
921 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
922END INTERFACE
923
926 MODULE PROCEDURE georef_coord_dist
927END INTERFACE
928
929#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
930#define ARRAYOF_TYPE arrayof_georef_coord_array
931!define ARRAYOF_ORIGEQ 0
932#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
933#include "arrayof_pre.F90"
934! from arrayof
936
937PRIVATE
939 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
940 georef_coord_array_polygon, georef_coord_array_multipoint, &
942#ifdef HAVE_SHAPELIB
944#endif
946 georef_coord_new, georef_coord_array_new
947
948CONTAINS
949
950#include "arrayof_post.F90"
951
952! ===================
953! == georef_coord ==
954! ===================
958FUNCTION georef_coord_new(x, y) RESULT(this)
959DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
960DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
961TYPE(georef_coord) :: this
962
965
966END FUNCTION georef_coord_new
967
968
969SUBROUTINE georef_coord_delete(this)
970TYPE(georef_coord),INTENT(inout) :: this
971
972this%x = dmiss
973this%y = dmiss
974
975END SUBROUTINE georef_coord_delete
976
977
978ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
979TYPE(georef_coord),INTENT(in) :: this
980LOGICAL :: res
981
982res = .NOT. this == georef_coord_miss
983
984END FUNCTION georef_coord_c_e
985
986
993ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
994TYPE(georef_coord),INTENT(in) :: this
995DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
996DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
997
998IF (PRESENT(x)) x = this%x
999IF (PRESENT(y)) y = this%y
1000
1001END SUBROUTINE georef_coord_getval
1002
1003
1012ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1013TYPE(georef_coord),INTENT(in) :: this
1014TYPE(geo_proj),INTENT(in) :: proj
1015DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1016DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1017DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1018DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1019
1020DOUBLE PRECISION :: llon, llat
1021
1022IF (PRESENT(x)) x = this%x
1023IF (PRESENT(y)) y = this%y
1024IF (PRESENT(lon) .OR. present(lat)) THEN
1026 IF (PRESENT(lon)) lon = llon
1027 IF (PRESENT(lat)) lat = llat
1028ENDIF
1029
1030END SUBROUTINE georef_coord_proj_getval
1031
1032
1033! document and improve
1034ELEMENTAL FUNCTION getlat(this)
1035TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1036DOUBLE PRECISION :: getlat ! latitudine geografica
1037
1038getlat = this%y ! change!!!
1039
1040END FUNCTION getlat
1041
1042! document and improve
1043ELEMENTAL FUNCTION getlon(this)
1044TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1045DOUBLE PRECISION :: getlon ! longitudine geografica
1046
1047getlon = this%x ! change!!!
1048
1049END FUNCTION getlon
1050
1051
1052ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1053TYPE(georef_coord),INTENT(IN) :: this, that
1054LOGICAL :: res
1055
1056res = (this%x == that%x .AND. this%y == that%y)
1057
1058END FUNCTION georef_coord_eq
1059
1060
1061ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1062TYPE(georef_coord),INTENT(IN) :: this, that
1063LOGICAL :: res
1064
1065res = (this%x >= that%x .AND. this%y >= that%y)
1066
1067END FUNCTION georef_coord_ge
1068
1069
1070ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1071TYPE(georef_coord),INTENT(IN) :: this, that
1072LOGICAL :: res
1073
1074res = (this%x <= that%x .AND. this%y <= that%y)
1075
1076END FUNCTION georef_coord_le
1077
1078
1079ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1080TYPE(georef_coord),INTENT(IN) :: this, that
1081LOGICAL :: res
1082
1083res = .NOT.(this == that)
1084
1085END FUNCTION georef_coord_ne
1086
1087
1093SUBROUTINE georef_coord_read_unit(this, unit)
1094TYPE(georef_coord),INTENT(out) :: this
1095INTEGER, INTENT(in) :: unit
1096
1097CALL georef_coord_vect_read_unit((/this/), unit)
1098
1099END SUBROUTINE georef_coord_read_unit
1100
1101
1107SUBROUTINE georef_coord_vect_read_unit(this, unit)
1108TYPE(georef_coord) :: this(:)
1109INTEGER, INTENT(in) :: unit
1110
1111CHARACTER(len=40) :: form
1112INTEGER :: i
1113
1114INQUIRE(unit, form=form)
1115IF (form == 'FORMATTED') THEN
1116 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1117!TODO bug gfortran compiler !
1118!missing values are unredeable when formatted
1119ELSE
1120 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1121ENDIF
1122
1123END SUBROUTINE georef_coord_vect_read_unit
1124
1125
1130SUBROUTINE georef_coord_write_unit(this, unit)
1131TYPE(georef_coord),INTENT(in) :: this
1132INTEGER, INTENT(in) :: unit
1133
1134CALL georef_coord_vect_write_unit((/this/), unit)
1135
1136END SUBROUTINE georef_coord_write_unit
1137
1138
1143SUBROUTINE georef_coord_vect_write_unit(this, unit)
1144TYPE(georef_coord),INTENT(in) :: this(:)
1145INTEGER, INTENT(in) :: unit
1146
1147CHARACTER(len=40) :: form
1148INTEGER :: i
1149
1150INQUIRE(unit, form=form)
1151IF (form == 'FORMATTED') THEN
1152 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1153!TODO bug gfortran compiler !
1154!missing values are unredeable when formatted
1155ELSE
1156 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1157ENDIF
1158
1159END SUBROUTINE georef_coord_vect_write_unit
1160
1161
1164FUNCTION georef_coord_dist(this, that) RESULT(dist)
1166TYPE(georef_coord), INTENT (IN) :: this
1167TYPE(georef_coord), INTENT (IN) :: that
1168DOUBLE PRECISION :: dist
1169
1170DOUBLE PRECISION :: x,y
1171! Distanza approssimata, valida per piccoli angoli
1172
1173x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1174y = (this%y-that%y)
1175dist = sqrt(x**2 + y**2)*degrad*rearth
1176
1177END FUNCTION georef_coord_dist
1178
1179
1185FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1186TYPE(georef_coord),INTENT(IN) :: this
1187TYPE(georef_coord),INTENT(IN) :: coordmin
1188TYPE(georef_coord),INTENT(IN) :: coordmax
1189LOGICAL :: res
1190
1191res = (this >= coordmin .AND. this <= coordmax)
1192
1193END FUNCTION georef_coord_inside_rectang
1194
1195
1196! ========================
1197! == georef_coord_array ==
1198! ========================
1204FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1205DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1206DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1207INTEGER,INTENT(in),OPTIONAL :: topo
1208TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1209TYPE(georef_coord_array) :: this
1210
1211INTEGER :: lsize
1212
1213IF (PRESENT(x) .AND. PRESENT(y)) THEN
1214 lsize = min(SIZE(x), SIZE(y))
1215 ALLOCATE(this%coord(lsize))
1216 this%coord(1:lsize)%x = x(1:lsize)
1217 this%coord(1:lsize)%y = y(1:lsize)
1218ENDIF
1219this%topo = optio_l(topo)
1221
1222END FUNCTION georef_coord_array_new
1223
1224
1225SUBROUTINE georef_coord_array_delete(this)
1226TYPE(georef_coord_array),INTENT(inout) :: this
1227
1228TYPE(georef_coord_array) :: lobj
1229
1230this = lobj
1231
1232END SUBROUTINE georef_coord_array_delete
1233
1234
1235ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1236TYPE(georef_coord_array),INTENT(in) :: this
1237LOGICAL :: res
1238
1239res = ALLOCATED(this%coord)
1240
1241END FUNCTION georef_coord_array_c_e
1242
1243
1248SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1249TYPE(georef_coord_array),INTENT(in) :: this
1250DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1251DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1252! allocatable per vedere di nascosto l'effetto che fa
1253INTEGER,OPTIONAL,INTENT(out) :: topo
1254TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1255
1256
1257IF (PRESENT(x)) THEN
1258 IF (ALLOCATED(this%coord)) THEN
1259 x = this%coord%x
1260 ENDIF
1261ENDIF
1262IF (PRESENT(y)) THEN
1263 IF (ALLOCATED(this%coord)) THEN
1264 y = this%coord%y
1265 ENDIF
1266ENDIF
1267IF (PRESENT(topo)) topo = this%topo
1269
1270END SUBROUTINE georef_coord_array_getval
1271
1272
1278SUBROUTINE georef_coord_array_compute_bbox(this)
1279TYPE(georef_coord_array),INTENT(inout) :: this
1280
1281IF (ALLOCATED(this%coord)) THEN
1282 this%bbox(1)%x = minval(this%coord(:)%x)
1283 this%bbox(1)%y = minval(this%coord(:)%y)
1284 this%bbox(2)%x = maxval(this%coord(:)%x)
1285 this%bbox(2)%y = maxval(this%coord(:)%y)
1286 this%bbox_updated = .true.
1287ENDIF
1288
1289END SUBROUTINE georef_coord_array_compute_bbox
1290
1291#ifdef HAVE_SHAPELIB
1292! internal method for importing a single shape
1293SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1294TYPE(georef_coord_array),INTENT(OUT) :: this
1295TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1296INTEGER,INTENT(IN) :: nshp
1297
1298TYPE(shpobject) :: shpobj
1299
1300! read shape object
1301shpobj = shpreadobject(shphandle, nshp)
1302IF (.NOT.shpisnull(shpobj)) THEN
1303! import it in georef_coord object
1304 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1305 topo=shpobj%nshptype)
1306 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1307 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1308 ELSE IF (ALLOCATED(this%parts)) THEN
1309 DEALLOCATE(this%parts)
1310 ENDIF
1311 CALL shpdestroyobject(shpobj)
1312 CALL compute_bbox(this)
1313ENDIF
1314
1315
1316END SUBROUTINE georef_coord_array_import
1317
1318
1319! internal method for exporting a single shape
1320SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1321TYPE(georef_coord_array),INTENT(in) :: this
1322TYPE(shpfileobject),INTENT(inout) :: shphandle
1323INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1324
1325INTEGER :: i
1326TYPE(shpobject) :: shpobj
1327
1328IF (ALLOCATED(this%coord)) THEN
1329 IF (ALLOCATED(this%parts)) THEN
1330 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1331 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1332 ELSE
1333 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1334 this%coord(:)%x, this%coord(:)%y)
1335 ENDIF
1336ELSE
1337 RETURN
1338ENDIF
1339
1340IF (.NOT.shpisnull(shpobj)) THEN
1341 i = shpwriteobject(shphandle, nshp, shpobj)
1342 CALL shpdestroyobject(shpobj)
1343ENDIF
1344
1345END SUBROUTINE georef_coord_array_export
1346
1347
1358SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1359TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1360CHARACTER(len=*),INTENT(in) :: shpfile
1361
1362REAL(kind=fp_d) :: minb(4), maxb(4)
1363INTEGER :: i, ns, shptype, dbfnf, dbfnr
1364TYPE(shpfileobject) :: shphandle
1365
1366shphandle = shpopen(trim(shpfile), 'rb')
1367IF (shpfileisnull(shphandle)) THEN
1368 ! log here
1369 CALL raise_error()
1370 RETURN
1371ENDIF
1372
1373! get info about file
1374CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1375IF (ns > 0) THEN ! allocate and read the object
1377 DO i = 1, ns
1378 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1379 ENDDO
1380ENDIF
1381
1382CALL shpclose(shphandle)
1383! pack object to save memory
1385
1386END SUBROUTINE arrayof_georef_coord_array_import
1387
1388
1394SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1395TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1396CHARACTER(len=*),INTENT(in) :: shpfile
1397
1398INTEGER :: i
1399TYPE(shpfileobject) :: shphandle
1400
1401IF (this%arraysize > 0) THEN
1402 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1403ELSE
1404 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1405ENDIF
1406IF (shpfileisnull(shphandle)) THEN
1407 ! log here
1408 CALL raise_error()
1409 RETURN
1410ENDIF
1411
1412DO i = 1, this%arraysize
1413 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1414ENDDO
1415
1416CALL shpclose(shphandle)
1417
1418END SUBROUTINE arrayof_georef_coord_array_export
1419#endif
1420
1432FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1433TYPE(georef_coord), INTENT(IN) :: this
1434TYPE(georef_coord_array), INTENT(IN) :: poly
1435LOGICAL :: inside
1436
1437INTEGER :: i
1438
1439inside = .false.
1441IF (.NOT.ALLOCATED(poly%coord)) RETURN
1442! if outside bounding box stop here
1443IF (poly%bbox_updated) THEN
1444 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1445ENDIF
1446
1447IF (ALLOCATED(poly%parts)) THEN
1448 DO i = 1, SIZE(poly%parts)-1
1450 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1451 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1452 ENDDO
1453 IF (SIZE(poly%parts) > 0) THEN ! safety check
1455 poly%coord(poly%parts(i)+1:)%x, &
1456 poly%coord(poly%parts(i)+1:)%y)
1457 ENDIF
1458
1459ELSE
1460 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1461 inside = pointinpoly(this%x, this%y, &
1462 poly%coord(:)%x, poly%coord(:)%y)
1463ENDIF
1464
1465CONTAINS
1466
1467FUNCTION pointinpoly(x, y, px, py)
1468DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1469LOGICAL :: pointinpoly
1470
1471INTEGER :: i, j, starti
1472
1473pointinpoly = .false.
1474
1475IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1476 starti = 2
1477 j = 1
1478ELSE ! unclosed polygon
1479 starti = 1
1480 j = SIZE(px)
1481ENDIF
1482DO i = starti, SIZE(px)
1483 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1484 (py(j) <= y .AND. y < py(i))) THEN
1485 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1486 pointinpoly = .NOT. pointinpoly
1487 ENDIF
1488 ENDIF
1489 j = i
1490ENDDO
1491
1492END FUNCTION pointinpoly
1493
1494END FUNCTION georef_coord_inside
1495
1496
1497
Compute forward coordinate transformation from geographical system to projected system. Definition: geo_proj_class.F90:295 Compute backward coordinate transformation from projected system to geographical system. Definition: geo_proj_class.F90:301 Quick method to append an element to the array. Definition: georef_coord_class.F90:401 Compute the distance in m between two points. Definition: georef_coord_class.F90:344 Export an array of georef_coord_array objects to a file in ESRI/Shapefile format. Definition: georef_coord_class.F90:321 Methods for returning the value of object members. Definition: georef_coord_class.F90:282 Import an array of georef_coord_array objects from a file in ESRI/Shapefile format. Definition: georef_coord_class.F90:315 Method for inserting elements of the array at a desired position. Definition: georef_coord_class.F90:392 Determine whether a point lies inside a polygon or a rectangle. Definition: georef_coord_class.F90:339 Method for packing the array object reducing at a minimum the memory occupation, without destroying i... Definition: georef_coord_class.F90:424 Read a single georef_coord object or an array of georef_coord objects from a Fortran FORMATTED or UNF... Definition: georef_coord_class.F90:328 Method for removing elements of the array at a desired position. Definition: georef_coord_class.F90:407 Write a single georef_coord object or an array of georef_coord objects to a Fortran FORMATTED or UNFO... Definition: georef_coord_class.F90:334 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:227 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:253 Derive type defining a single georeferenced point, either in geodetic or in projected coordinates. Definition: georef_coord_class.F90:241 |