libsim Versione 7.2.4
|
◆ georef_coord_read_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.
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
887USE geo_proj_class
888#ifdef HAVE_SHAPELIB
889USE shapelib
890#endif
891IMPLICIT NONE
892
898 PRIVATE
899 DOUBLE PRECISION :: x=dmiss, y=dmiss
901
904
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.
918
919INTEGER,PARAMETER :: georef_coord_array_point = 1
920INTEGER,PARAMETER :: georef_coord_array_arc = 3
921INTEGER,PARAMETER :: georef_coord_array_polygon = 5
922INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
923
924
929 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
930END INTERFACE
931
934 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
935END INTERFACE
936
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
947INTERFACE OPERATOR (==)
948 MODULE PROCEDURE georef_coord_eq
949END INTERFACE
950
952INTERFACE OPERATOR (/=)
953 MODULE PROCEDURE georef_coord_ne
954END INTERFACE
955
958INTERFACE OPERATOR (>=)
959 MODULE PROCEDURE georef_coord_ge
960END INTERFACE
961
964INTERFACE OPERATOR (<=)
965 MODULE PROCEDURE georef_coord_le
966END INTERFACE
967
968#ifdef HAVE_SHAPELIB
969
971INTERFACE import
972 MODULE PROCEDURE arrayof_georef_coord_array_import
973END INTERFACE
974
978 MODULE PROCEDURE arrayof_georef_coord_array_export
979END INTERFACE
980#endif
981
985 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
986END INTERFACE
987
991 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
992END INTERFACE
993
996 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
997END INTERFACE
998
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
1011
1012PRIVATE
1014 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1015 georef_coord_array_polygon, georef_coord_array_multipoint, &
1017#ifdef HAVE_SHAPELIB
1019#endif
1021 georef_coord_new, georef_coord_array_new
1022
1023CONTAINS
1024
1025#include "arrayof_post.F90"
1026
1027! ===================
1028! == georef_coord ==
1029! ===================
1033FUNCTION georef_coord_new(x, y) RESULT(this)
1034DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1035DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1036TYPE(georef_coord) :: this
1037
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
1068ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1069TYPE(georef_coord),INTENT(in) :: this
1070DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1071DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1072
1073IF (PRESENT(x)) x = this%x
1074IF (PRESENT(y)) y = this%y
1075
1076END SUBROUTINE georef_coord_getval
1077
1078
1087ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1088TYPE(georef_coord),INTENT(in) :: this
1089TYPE(geo_proj),INTENT(in) :: proj
1090DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1091DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1092DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1093DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
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
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
1168SUBROUTINE georef_coord_read_unit(this, unit)
1169TYPE(georef_coord),INTENT(out) :: this
1170INTEGER, INTENT(in) :: unit
1171
1172CALL georef_coord_vect_read_unit((/this/), unit)
1173
1174END SUBROUTINE georef_coord_read_unit
1175
1176
1182SUBROUTINE georef_coord_vect_read_unit(this, unit)
1183TYPE(georef_coord) :: this(:)
1184INTEGER, INTENT(in) :: unit
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
1205SUBROUTINE georef_coord_write_unit(this, unit)
1206TYPE(georef_coord),INTENT(in) :: this
1207INTEGER, INTENT(in) :: unit
1208
1209CALL georef_coord_vect_write_unit((/this/), unit)
1210
1211END SUBROUTINE georef_coord_write_unit
1212
1213
1218SUBROUTINE georef_coord_vect_write_unit(this, unit)
1219TYPE(georef_coord),INTENT(in) :: this(:)
1220INTEGER, INTENT(in) :: unit
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
1239FUNCTION georef_coord_dist(this, that) RESULT(dist)
1241TYPE(georef_coord), INTENT (IN) :: this
1242TYPE(georef_coord), INTENT (IN) :: that
1243DOUBLE PRECISION :: dist
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
1260FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1261TYPE(georef_coord),INTENT(IN) :: this
1262TYPE(georef_coord),INTENT(IN) :: coordmin
1263TYPE(georef_coord),INTENT(IN) :: coordmax
1264LOGICAL :: res
1265
1266res = (this >= coordmin .AND. this <= coordmax)
1267
1268END FUNCTION georef_coord_inside_rectang
1269
1270
1271! ========================
1272! == georef_coord_array ==
1273! ========================
1279FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1280DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1281DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1282INTEGER,INTENT(in),OPTIONAL :: topo
1283TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
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)
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
1323SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1324TYPE(georef_coord_array),INTENT(in) :: this
1325DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1326DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1327! allocatable per vedere di nascosto l'effetto che fa
1328INTEGER,OPTIONAL,INTENT(out) :: topo
1329TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
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
1344
1345END SUBROUTINE georef_coord_array_getval
1346
1347
1353SUBROUTINE georef_coord_array_compute_bbox(this)
1354TYPE(georef_coord_array),INTENT(inout) :: this
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
1433SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1434TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1435CHARACTER(len=*),INTENT(in) :: shpfile
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
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
1460
1461END SUBROUTINE arrayof_georef_coord_array_import
1462
1463
1469SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1470TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1471CHARACTER(len=*),INTENT(in) :: shpfile
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
1507FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1508TYPE(georef_coord), INTENT(IN) :: this
1509TYPE(georef_coord_array), INTENT(IN) :: poly
1510LOGICAL :: inside
1511
1512INTEGER :: i
1513
1514inside = .false.
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
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
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
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 |