libsim Versione 7.1.11
|
◆ 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 856 del file georef_coord_class.F90. 857! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
858! authors:
859! Davide Cesari <dcesari@arpa.emr.it>
860! Paolo Patruno <ppatruno@arpa.emr.it>
861
862! This program is free software; you can redistribute it and/or
863! modify it under the terms of the GNU General Public License as
864! published by the Free Software Foundation; either version 2 of
865! the License, or (at your option) any later version.
866
867! This program is distributed in the hope that it will be useful,
868! but WITHOUT ANY WARRANTY; without even the implied warranty of
869! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
870! GNU General Public License for more details.
871
872! You should have received a copy of the GNU General Public License
873! along with this program. If not, see <http://www.gnu.org/licenses/>.
874#include "config.h"
875
893USE geo_proj_class
894#ifdef HAVE_SHAPELIB
895USE shapelib
896#endif
897IMPLICIT NONE
898
904 PRIVATE
905 DOUBLE PRECISION :: x=dmiss, y=dmiss
907
910
916 PRIVATE
917 INTEGER,ALLOCATABLE :: parts(:)
918 TYPE(georef_coord),ALLOCATABLE :: coord(:)
919 INTEGER :: topo=imiss
920 TYPE(geo_proj) :: proj
921 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
922 LOGICAL :: bbox_updated=.false.
924
925INTEGER,PARAMETER :: georef_coord_array_point = 1
926INTEGER,PARAMETER :: georef_coord_array_arc = 3
927INTEGER,PARAMETER :: georef_coord_array_polygon = 5
928INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
929
930
935 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
936END INTERFACE
937
940 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
941END INTERFACE
942
945 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
946END INTERFACE
947
948INTERFACE compute_bbox
949 MODULE PROCEDURE georef_coord_array_compute_bbox
950END INTERFACE
951
953INTERFACE OPERATOR (==)
954 MODULE PROCEDURE georef_coord_eq
955END INTERFACE
956
958INTERFACE OPERATOR (/=)
959 MODULE PROCEDURE georef_coord_ne
960END INTERFACE
961
964INTERFACE OPERATOR (>=)
965 MODULE PROCEDURE georef_coord_ge
966END INTERFACE
967
970INTERFACE OPERATOR (<=)
971 MODULE PROCEDURE georef_coord_le
972END INTERFACE
973
974#ifdef HAVE_SHAPELIB
975
977INTERFACE import
978 MODULE PROCEDURE arrayof_georef_coord_array_import
979END INTERFACE
980
984 MODULE PROCEDURE arrayof_georef_coord_array_export
985END INTERFACE
986#endif
987
991 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
992END INTERFACE
993
997 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
998END INTERFACE
999
1002 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1003END INTERFACE
1004
1007 MODULE PROCEDURE georef_coord_dist
1008END INTERFACE
1009
1010#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1011#define ARRAYOF_TYPE arrayof_georef_coord_array
1012!define ARRAYOF_ORIGEQ 0
1013#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1014#include "arrayof_pre.F90"
1015! from arrayof
1017
1018PRIVATE
1020 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1021 georef_coord_array_polygon, georef_coord_array_multipoint, &
1023#ifdef HAVE_SHAPELIB
1025#endif
1027 georef_coord_new, georef_coord_array_new
1028
1029CONTAINS
1030
1031#include "arrayof_post.F90"
1032
1033! ===================
1034! == georef_coord ==
1035! ===================
1039FUNCTION georef_coord_new(x, y) RESULT(this)
1040DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1041DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1042TYPE(georef_coord) :: this
1043
1046
1047END FUNCTION georef_coord_new
1048
1049
1050SUBROUTINE georef_coord_delete(this)
1051TYPE(georef_coord),INTENT(inout) :: this
1052
1053this%x = dmiss
1054this%y = dmiss
1055
1056END SUBROUTINE georef_coord_delete
1057
1058
1059ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1060TYPE(georef_coord),INTENT(in) :: this
1061LOGICAL :: res
1062
1063res = .NOT. this == georef_coord_miss
1064
1065END FUNCTION georef_coord_c_e
1066
1067
1074ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1075TYPE(georef_coord),INTENT(in) :: this
1076DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1077DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1078
1079IF (PRESENT(x)) x = this%x
1080IF (PRESENT(y)) y = this%y
1081
1082END SUBROUTINE georef_coord_getval
1083
1084
1093ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1094TYPE(georef_coord),INTENT(in) :: this
1095TYPE(geo_proj),INTENT(in) :: proj
1096DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1097DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1098DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1099DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1100
1101DOUBLE PRECISION :: llon, llat
1102
1103IF (PRESENT(x)) x = this%x
1104IF (PRESENT(y)) y = this%y
1105IF (PRESENT(lon) .OR. present(lat)) THEN
1107 IF (PRESENT(lon)) lon = llon
1108 IF (PRESENT(lat)) lat = llat
1109ENDIF
1110
1111END SUBROUTINE georef_coord_proj_getval
1112
1113
1114! document and improve
1115ELEMENTAL FUNCTION getlat(this)
1116TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1117DOUBLE PRECISION :: getlat ! latitudine geografica
1118
1119getlat = this%y ! change!!!
1120
1121END FUNCTION getlat
1122
1123! document and improve
1124ELEMENTAL FUNCTION getlon(this)
1125TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1126DOUBLE PRECISION :: getlon ! longitudine geografica
1127
1128getlon = this%x ! change!!!
1129
1130END FUNCTION getlon
1131
1132
1133ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1134TYPE(georef_coord),INTENT(IN) :: this, that
1135LOGICAL :: res
1136
1137res = (this%x == that%x .AND. this%y == that%y)
1138
1139END FUNCTION georef_coord_eq
1140
1141
1142ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1143TYPE(georef_coord),INTENT(IN) :: this, that
1144LOGICAL :: res
1145
1146res = (this%x >= that%x .AND. this%y >= that%y)
1147
1148END FUNCTION georef_coord_ge
1149
1150
1151ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1152TYPE(georef_coord),INTENT(IN) :: this, that
1153LOGICAL :: res
1154
1155res = (this%x <= that%x .AND. this%y <= that%y)
1156
1157END FUNCTION georef_coord_le
1158
1159
1160ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1161TYPE(georef_coord),INTENT(IN) :: this, that
1162LOGICAL :: res
1163
1164res = .NOT.(this == that)
1165
1166END FUNCTION georef_coord_ne
1167
1168
1174SUBROUTINE georef_coord_read_unit(this, unit)
1175TYPE(georef_coord),INTENT(out) :: this
1176INTEGER, INTENT(in) :: unit
1177
1178CALL georef_coord_vect_read_unit((/this/), unit)
1179
1180END SUBROUTINE georef_coord_read_unit
1181
1182
1188SUBROUTINE georef_coord_vect_read_unit(this, unit)
1189TYPE(georef_coord) :: this(:)
1190INTEGER, INTENT(in) :: unit
1191
1192CHARACTER(len=40) :: form
1193INTEGER :: i
1194
1195INQUIRE(unit, form=form)
1196IF (form == 'FORMATTED') THEN
1197 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1198!TODO bug gfortran compiler !
1199!missing values are unredeable when formatted
1200ELSE
1201 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1202ENDIF
1203
1204END SUBROUTINE georef_coord_vect_read_unit
1205
1206
1211SUBROUTINE georef_coord_write_unit(this, unit)
1212TYPE(georef_coord),INTENT(in) :: this
1213INTEGER, INTENT(in) :: unit
1214
1215CALL georef_coord_vect_write_unit((/this/), unit)
1216
1217END SUBROUTINE georef_coord_write_unit
1218
1219
1224SUBROUTINE georef_coord_vect_write_unit(this, unit)
1225TYPE(georef_coord),INTENT(in) :: this(:)
1226INTEGER, INTENT(in) :: unit
1227
1228CHARACTER(len=40) :: form
1229INTEGER :: i
1230
1231INQUIRE(unit, form=form)
1232IF (form == 'FORMATTED') THEN
1233 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1234!TODO bug gfortran compiler !
1235!missing values are unredeable when formatted
1236ELSE
1237 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1238ENDIF
1239
1240END SUBROUTINE georef_coord_vect_write_unit
1241
1242
1245FUNCTION georef_coord_dist(this, that) RESULT(dist)
1247TYPE(georef_coord), INTENT (IN) :: this
1248TYPE(georef_coord), INTENT (IN) :: that
1249DOUBLE PRECISION :: dist
1250
1251DOUBLE PRECISION :: x,y
1252! Distanza approssimata, valida per piccoli angoli
1253
1254x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1255y = (this%y-that%y)
1256dist = sqrt(x**2 + y**2)*degrad*rearth
1257
1258END FUNCTION georef_coord_dist
1259
1260
1266FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1267TYPE(georef_coord),INTENT(IN) :: this
1268TYPE(georef_coord),INTENT(IN) :: coordmin
1269TYPE(georef_coord),INTENT(IN) :: coordmax
1270LOGICAL :: res
1271
1272res = (this >= coordmin .AND. this <= coordmax)
1273
1274END FUNCTION georef_coord_inside_rectang
1275
1276
1277! ========================
1278! == georef_coord_array ==
1279! ========================
1285FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1286DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1287DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1288INTEGER,INTENT(in),OPTIONAL :: topo
1289TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1290TYPE(georef_coord_array) :: this
1291
1292INTEGER :: lsize
1293
1294IF (PRESENT(x) .AND. PRESENT(y)) THEN
1295 lsize = min(SIZE(x), SIZE(y))
1296 ALLOCATE(this%coord(lsize))
1297 this%coord(1:lsize)%x = x(1:lsize)
1298 this%coord(1:lsize)%y = y(1:lsize)
1299ENDIF
1300this%topo = optio_l(topo)
1302
1303END FUNCTION georef_coord_array_new
1304
1305
1306SUBROUTINE georef_coord_array_delete(this)
1307TYPE(georef_coord_array),INTENT(inout) :: this
1308
1309TYPE(georef_coord_array) :: lobj
1310
1311this = lobj
1312
1313END SUBROUTINE georef_coord_array_delete
1314
1315
1316ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1317TYPE(georef_coord_array),INTENT(in) :: this
1318LOGICAL :: res
1319
1320res = ALLOCATED(this%coord)
1321
1322END FUNCTION georef_coord_array_c_e
1323
1324
1329SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1330TYPE(georef_coord_array),INTENT(in) :: this
1331DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1332DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1333! allocatable per vedere di nascosto l'effetto che fa
1334INTEGER,OPTIONAL,INTENT(out) :: topo
1335TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1336
1337
1338IF (PRESENT(x)) THEN
1339 IF (ALLOCATED(this%coord)) THEN
1340 x = this%coord%x
1341 ENDIF
1342ENDIF
1343IF (PRESENT(y)) THEN
1344 IF (ALLOCATED(this%coord)) THEN
1345 y = this%coord%y
1346 ENDIF
1347ENDIF
1348IF (PRESENT(topo)) topo = this%topo
1350
1351END SUBROUTINE georef_coord_array_getval
1352
1353
1359SUBROUTINE georef_coord_array_compute_bbox(this)
1360TYPE(georef_coord_array),INTENT(inout) :: this
1361
1362IF (ALLOCATED(this%coord)) THEN
1363 this%bbox(1)%x = minval(this%coord(:)%x)
1364 this%bbox(1)%y = minval(this%coord(:)%y)
1365 this%bbox(2)%x = maxval(this%coord(:)%x)
1366 this%bbox(2)%y = maxval(this%coord(:)%y)
1367 this%bbox_updated = .true.
1368ENDIF
1369
1370END SUBROUTINE georef_coord_array_compute_bbox
1371
1372#ifdef HAVE_SHAPELIB
1373! internal method for importing a single shape
1374SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1375TYPE(georef_coord_array),INTENT(OUT) :: this
1376TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1377INTEGER,INTENT(IN) :: nshp
1378
1379TYPE(shpobject) :: shpobj
1380
1381! read shape object
1382shpobj = shpreadobject(shphandle, nshp)
1383IF (.NOT.shpisnull(shpobj)) THEN
1384! import it in georef_coord object
1385 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1386 topo=shpobj%nshptype)
1387 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1388 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1389 ELSE IF (ALLOCATED(this%parts)) THEN
1390 DEALLOCATE(this%parts)
1391 ENDIF
1392 CALL shpdestroyobject(shpobj)
1393 CALL compute_bbox(this)
1394ENDIF
1395
1396
1397END SUBROUTINE georef_coord_array_import
1398
1399
1400! internal method for exporting a single shape
1401SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1402TYPE(georef_coord_array),INTENT(in) :: this
1403TYPE(shpfileobject),INTENT(inout) :: shphandle
1404INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1405
1406INTEGER :: i
1407TYPE(shpobject) :: shpobj
1408
1409IF (ALLOCATED(this%coord)) THEN
1410 IF (ALLOCATED(this%parts)) THEN
1411 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1412 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1413 ELSE
1414 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1415 this%coord(:)%x, this%coord(:)%y)
1416 ENDIF
1417ELSE
1418 RETURN
1419ENDIF
1420
1421IF (.NOT.shpisnull(shpobj)) THEN
1422 i = shpwriteobject(shphandle, nshp, shpobj)
1423 CALL shpdestroyobject(shpobj)
1424ENDIF
1425
1426END SUBROUTINE georef_coord_array_export
1427
1428
1439SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1440TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1441CHARACTER(len=*),INTENT(in) :: shpfile
1442
1443REAL(kind=fp_d) :: minb(4), maxb(4)
1444INTEGER :: i, ns, shptype, dbfnf, dbfnr
1445TYPE(shpfileobject) :: shphandle
1446
1447shphandle = shpopen(trim(shpfile), 'rb')
1448IF (shpfileisnull(shphandle)) THEN
1449 ! log here
1450 CALL raise_error()
1451 RETURN
1452ENDIF
1453
1454! get info about file
1455CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1456IF (ns > 0) THEN ! allocate and read the object
1458 DO i = 1, ns
1459 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1460 ENDDO
1461ENDIF
1462
1463CALL shpclose(shphandle)
1464! pack object to save memory
1466
1467END SUBROUTINE arrayof_georef_coord_array_import
1468
1469
1475SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1476TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1477CHARACTER(len=*),INTENT(in) :: shpfile
1478
1479INTEGER :: i
1480TYPE(shpfileobject) :: shphandle
1481
1482IF (this%arraysize > 0) THEN
1483 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1484ELSE
1485 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1486ENDIF
1487IF (shpfileisnull(shphandle)) THEN
1488 ! log here
1489 CALL raise_error()
1490 RETURN
1491ENDIF
1492
1493DO i = 1, this%arraysize
1494 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1495ENDDO
1496
1497CALL shpclose(shphandle)
1498
1499END SUBROUTINE arrayof_georef_coord_array_export
1500#endif
1501
1513FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1514TYPE(georef_coord), INTENT(IN) :: this
1515TYPE(georef_coord_array), INTENT(IN) :: poly
1516LOGICAL :: inside
1517
1518INTEGER :: i
1519
1520inside = .false.
1522IF (.NOT.ALLOCATED(poly%coord)) RETURN
1523! if outside bounding box stop here
1524IF (poly%bbox_updated) THEN
1525 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1526ENDIF
1527
1528IF (ALLOCATED(poly%parts)) THEN
1529 DO i = 1, SIZE(poly%parts)-1
1531 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1532 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1533 ENDDO
1534 IF (SIZE(poly%parts) > 0) THEN ! safety check
1536 poly%coord(poly%parts(i)+1:)%x, &
1537 poly%coord(poly%parts(i)+1:)%y)
1538 ENDIF
1539
1540ELSE
1541 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1542 inside = pointinpoly(this%x, this%y, &
1543 poly%coord(:)%x, poly%coord(:)%y)
1544ENDIF
1545
1546CONTAINS
1547
1548FUNCTION pointinpoly(x, y, px, py)
1549DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1550LOGICAL :: pointinpoly
1551
1552INTEGER :: i, j, starti
1553
1554pointinpoly = .false.
1555
1556IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1557 starti = 2
1558 j = 1
1559ELSE ! unclosed polygon
1560 starti = 1
1561 j = SIZE(px)
1562ENDIF
1563DO i = starti, SIZE(px)
1564 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1565 (py(j) <= y .AND. y < py(i))) THEN
1566 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1567 pointinpoly = .NOT. pointinpoly
1568 ENDIF
1569 ENDIF
1570 j = i
1571ENDDO
1572
1573END FUNCTION pointinpoly
1574
1575END FUNCTION georef_coord_inside
1576
1577
1578
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 |