libsim Versione 7.1.11
|
◆ georef_coord_vect_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 870 del file georef_coord_class.F90. 871! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
872! authors:
873! Davide Cesari <dcesari@arpa.emr.it>
874! Paolo Patruno <ppatruno@arpa.emr.it>
875
876! This program is free software; you can redistribute it and/or
877! modify it under the terms of the GNU General Public License as
878! published by the Free Software Foundation; either version 2 of
879! the License, or (at your option) any later version.
880
881! This program is distributed in the hope that it will be useful,
882! but WITHOUT ANY WARRANTY; without even the implied warranty of
883! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
884! GNU General Public License for more details.
885
886! You should have received a copy of the GNU General Public License
887! along with this program. If not, see <http://www.gnu.org/licenses/>.
888#include "config.h"
889
907USE geo_proj_class
908#ifdef HAVE_SHAPELIB
909USE shapelib
910#endif
911IMPLICIT NONE
912
918 PRIVATE
919 DOUBLE PRECISION :: x=dmiss, y=dmiss
921
924
930 PRIVATE
931 INTEGER,ALLOCATABLE :: parts(:)
932 TYPE(georef_coord),ALLOCATABLE :: coord(:)
933 INTEGER :: topo=imiss
934 TYPE(geo_proj) :: proj
935 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
936 LOGICAL :: bbox_updated=.false.
938
939INTEGER,PARAMETER :: georef_coord_array_point = 1
940INTEGER,PARAMETER :: georef_coord_array_arc = 3
941INTEGER,PARAMETER :: georef_coord_array_polygon = 5
942INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
943
944
949 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
950END INTERFACE
951
954 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
955END INTERFACE
956
959 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
960END INTERFACE
961
962INTERFACE compute_bbox
963 MODULE PROCEDURE georef_coord_array_compute_bbox
964END INTERFACE
965
967INTERFACE OPERATOR (==)
968 MODULE PROCEDURE georef_coord_eq
969END INTERFACE
970
972INTERFACE OPERATOR (/=)
973 MODULE PROCEDURE georef_coord_ne
974END INTERFACE
975
978INTERFACE OPERATOR (>=)
979 MODULE PROCEDURE georef_coord_ge
980END INTERFACE
981
984INTERFACE OPERATOR (<=)
985 MODULE PROCEDURE georef_coord_le
986END INTERFACE
987
988#ifdef HAVE_SHAPELIB
989
991INTERFACE import
992 MODULE PROCEDURE arrayof_georef_coord_array_import
993END INTERFACE
994
998 MODULE PROCEDURE arrayof_georef_coord_array_export
999END INTERFACE
1000#endif
1001
1005 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1006END INTERFACE
1007
1011 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1012END INTERFACE
1013
1016 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1017END INTERFACE
1018
1021 MODULE PROCEDURE georef_coord_dist
1022END INTERFACE
1023
1024#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1025#define ARRAYOF_TYPE arrayof_georef_coord_array
1026!define ARRAYOF_ORIGEQ 0
1027#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1028#include "arrayof_pre.F90"
1029! from arrayof
1031
1032PRIVATE
1034 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1035 georef_coord_array_polygon, georef_coord_array_multipoint, &
1037#ifdef HAVE_SHAPELIB
1039#endif
1041 georef_coord_new, georef_coord_array_new
1042
1043CONTAINS
1044
1045#include "arrayof_post.F90"
1046
1047! ===================
1048! == georef_coord ==
1049! ===================
1053FUNCTION georef_coord_new(x, y) RESULT(this)
1054DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1055DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1056TYPE(georef_coord) :: this
1057
1060
1061END FUNCTION georef_coord_new
1062
1063
1064SUBROUTINE georef_coord_delete(this)
1065TYPE(georef_coord),INTENT(inout) :: this
1066
1067this%x = dmiss
1068this%y = dmiss
1069
1070END SUBROUTINE georef_coord_delete
1071
1072
1073ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1074TYPE(georef_coord),INTENT(in) :: this
1075LOGICAL :: res
1076
1077res = .NOT. this == georef_coord_miss
1078
1079END FUNCTION georef_coord_c_e
1080
1081
1088ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1089TYPE(georef_coord),INTENT(in) :: this
1090DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1091DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1092
1093IF (PRESENT(x)) x = this%x
1094IF (PRESENT(y)) y = this%y
1095
1096END SUBROUTINE georef_coord_getval
1097
1098
1107ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1108TYPE(georef_coord),INTENT(in) :: this
1109TYPE(geo_proj),INTENT(in) :: proj
1110DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1111DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1112DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1113DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1114
1115DOUBLE PRECISION :: llon, llat
1116
1117IF (PRESENT(x)) x = this%x
1118IF (PRESENT(y)) y = this%y
1119IF (PRESENT(lon) .OR. present(lat)) THEN
1121 IF (PRESENT(lon)) lon = llon
1122 IF (PRESENT(lat)) lat = llat
1123ENDIF
1124
1125END SUBROUTINE georef_coord_proj_getval
1126
1127
1128! document and improve
1129ELEMENTAL FUNCTION getlat(this)
1130TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1131DOUBLE PRECISION :: getlat ! latitudine geografica
1132
1133getlat = this%y ! change!!!
1134
1135END FUNCTION getlat
1136
1137! document and improve
1138ELEMENTAL FUNCTION getlon(this)
1139TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1140DOUBLE PRECISION :: getlon ! longitudine geografica
1141
1142getlon = this%x ! change!!!
1143
1144END FUNCTION getlon
1145
1146
1147ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1148TYPE(georef_coord),INTENT(IN) :: this, that
1149LOGICAL :: res
1150
1151res = (this%x == that%x .AND. this%y == that%y)
1152
1153END FUNCTION georef_coord_eq
1154
1155
1156ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1157TYPE(georef_coord),INTENT(IN) :: this, that
1158LOGICAL :: res
1159
1160res = (this%x >= that%x .AND. this%y >= that%y)
1161
1162END FUNCTION georef_coord_ge
1163
1164
1165ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1166TYPE(georef_coord),INTENT(IN) :: this, that
1167LOGICAL :: res
1168
1169res = (this%x <= that%x .AND. this%y <= that%y)
1170
1171END FUNCTION georef_coord_le
1172
1173
1174ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1175TYPE(georef_coord),INTENT(IN) :: this, that
1176LOGICAL :: res
1177
1178res = .NOT.(this == that)
1179
1180END FUNCTION georef_coord_ne
1181
1182
1188SUBROUTINE georef_coord_read_unit(this, unit)
1189TYPE(georef_coord),INTENT(out) :: this
1190INTEGER, INTENT(in) :: unit
1191
1192CALL georef_coord_vect_read_unit((/this/), unit)
1193
1194END SUBROUTINE georef_coord_read_unit
1195
1196
1202SUBROUTINE georef_coord_vect_read_unit(this, unit)
1203TYPE(georef_coord) :: this(:)
1204INTEGER, INTENT(in) :: unit
1205
1206CHARACTER(len=40) :: form
1207INTEGER :: i
1208
1209INQUIRE(unit, form=form)
1210IF (form == 'FORMATTED') THEN
1211 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1212!TODO bug gfortran compiler !
1213!missing values are unredeable when formatted
1214ELSE
1215 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1216ENDIF
1217
1218END SUBROUTINE georef_coord_vect_read_unit
1219
1220
1225SUBROUTINE georef_coord_write_unit(this, unit)
1226TYPE(georef_coord),INTENT(in) :: this
1227INTEGER, INTENT(in) :: unit
1228
1229CALL georef_coord_vect_write_unit((/this/), unit)
1230
1231END SUBROUTINE georef_coord_write_unit
1232
1233
1238SUBROUTINE georef_coord_vect_write_unit(this, unit)
1239TYPE(georef_coord),INTENT(in) :: this(:)
1240INTEGER, INTENT(in) :: unit
1241
1242CHARACTER(len=40) :: form
1243INTEGER :: i
1244
1245INQUIRE(unit, form=form)
1246IF (form == 'FORMATTED') THEN
1247 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1248!TODO bug gfortran compiler !
1249!missing values are unredeable when formatted
1250ELSE
1251 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1252ENDIF
1253
1254END SUBROUTINE georef_coord_vect_write_unit
1255
1256
1259FUNCTION georef_coord_dist(this, that) RESULT(dist)
1261TYPE(georef_coord), INTENT (IN) :: this
1262TYPE(georef_coord), INTENT (IN) :: that
1263DOUBLE PRECISION :: dist
1264
1265DOUBLE PRECISION :: x,y
1266! Distanza approssimata, valida per piccoli angoli
1267
1268x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1269y = (this%y-that%y)
1270dist = sqrt(x**2 + y**2)*degrad*rearth
1271
1272END FUNCTION georef_coord_dist
1273
1274
1280FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1281TYPE(georef_coord),INTENT(IN) :: this
1282TYPE(georef_coord),INTENT(IN) :: coordmin
1283TYPE(georef_coord),INTENT(IN) :: coordmax
1284LOGICAL :: res
1285
1286res = (this >= coordmin .AND. this <= coordmax)
1287
1288END FUNCTION georef_coord_inside_rectang
1289
1290
1291! ========================
1292! == georef_coord_array ==
1293! ========================
1299FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1300DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1301DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1302INTEGER,INTENT(in),OPTIONAL :: topo
1303TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1304TYPE(georef_coord_array) :: this
1305
1306INTEGER :: lsize
1307
1308IF (PRESENT(x) .AND. PRESENT(y)) THEN
1309 lsize = min(SIZE(x), SIZE(y))
1310 ALLOCATE(this%coord(lsize))
1311 this%coord(1:lsize)%x = x(1:lsize)
1312 this%coord(1:lsize)%y = y(1:lsize)
1313ENDIF
1314this%topo = optio_l(topo)
1316
1317END FUNCTION georef_coord_array_new
1318
1319
1320SUBROUTINE georef_coord_array_delete(this)
1321TYPE(georef_coord_array),INTENT(inout) :: this
1322
1323TYPE(georef_coord_array) :: lobj
1324
1325this = lobj
1326
1327END SUBROUTINE georef_coord_array_delete
1328
1329
1330ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1331TYPE(georef_coord_array),INTENT(in) :: this
1332LOGICAL :: res
1333
1334res = ALLOCATED(this%coord)
1335
1336END FUNCTION georef_coord_array_c_e
1337
1338
1343SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1344TYPE(georef_coord_array),INTENT(in) :: this
1345DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1346DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1347! allocatable per vedere di nascosto l'effetto che fa
1348INTEGER,OPTIONAL,INTENT(out) :: topo
1349TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1350
1351
1352IF (PRESENT(x)) THEN
1353 IF (ALLOCATED(this%coord)) THEN
1354 x = this%coord%x
1355 ENDIF
1356ENDIF
1357IF (PRESENT(y)) THEN
1358 IF (ALLOCATED(this%coord)) THEN
1359 y = this%coord%y
1360 ENDIF
1361ENDIF
1362IF (PRESENT(topo)) topo = this%topo
1364
1365END SUBROUTINE georef_coord_array_getval
1366
1367
1373SUBROUTINE georef_coord_array_compute_bbox(this)
1374TYPE(georef_coord_array),INTENT(inout) :: this
1375
1376IF (ALLOCATED(this%coord)) THEN
1377 this%bbox(1)%x = minval(this%coord(:)%x)
1378 this%bbox(1)%y = minval(this%coord(:)%y)
1379 this%bbox(2)%x = maxval(this%coord(:)%x)
1380 this%bbox(2)%y = maxval(this%coord(:)%y)
1381 this%bbox_updated = .true.
1382ENDIF
1383
1384END SUBROUTINE georef_coord_array_compute_bbox
1385
1386#ifdef HAVE_SHAPELIB
1387! internal method for importing a single shape
1388SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1389TYPE(georef_coord_array),INTENT(OUT) :: this
1390TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1391INTEGER,INTENT(IN) :: nshp
1392
1393TYPE(shpobject) :: shpobj
1394
1395! read shape object
1396shpobj = shpreadobject(shphandle, nshp)
1397IF (.NOT.shpisnull(shpobj)) THEN
1398! import it in georef_coord object
1399 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1400 topo=shpobj%nshptype)
1401 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1402 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1403 ELSE IF (ALLOCATED(this%parts)) THEN
1404 DEALLOCATE(this%parts)
1405 ENDIF
1406 CALL shpdestroyobject(shpobj)
1407 CALL compute_bbox(this)
1408ENDIF
1409
1410
1411END SUBROUTINE georef_coord_array_import
1412
1413
1414! internal method for exporting a single shape
1415SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1416TYPE(georef_coord_array),INTENT(in) :: this
1417TYPE(shpfileobject),INTENT(inout) :: shphandle
1418INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1419
1420INTEGER :: i
1421TYPE(shpobject) :: shpobj
1422
1423IF (ALLOCATED(this%coord)) THEN
1424 IF (ALLOCATED(this%parts)) THEN
1425 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1426 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1427 ELSE
1428 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1429 this%coord(:)%x, this%coord(:)%y)
1430 ENDIF
1431ELSE
1432 RETURN
1433ENDIF
1434
1435IF (.NOT.shpisnull(shpobj)) THEN
1436 i = shpwriteobject(shphandle, nshp, shpobj)
1437 CALL shpdestroyobject(shpobj)
1438ENDIF
1439
1440END SUBROUTINE georef_coord_array_export
1441
1442
1453SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1454TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1455CHARACTER(len=*),INTENT(in) :: shpfile
1456
1457REAL(kind=fp_d) :: minb(4), maxb(4)
1458INTEGER :: i, ns, shptype, dbfnf, dbfnr
1459TYPE(shpfileobject) :: shphandle
1460
1461shphandle = shpopen(trim(shpfile), 'rb')
1462IF (shpfileisnull(shphandle)) THEN
1463 ! log here
1464 CALL raise_error()
1465 RETURN
1466ENDIF
1467
1468! get info about file
1469CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1470IF (ns > 0) THEN ! allocate and read the object
1472 DO i = 1, ns
1473 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1474 ENDDO
1475ENDIF
1476
1477CALL shpclose(shphandle)
1478! pack object to save memory
1480
1481END SUBROUTINE arrayof_georef_coord_array_import
1482
1483
1489SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1490TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1491CHARACTER(len=*),INTENT(in) :: shpfile
1492
1493INTEGER :: i
1494TYPE(shpfileobject) :: shphandle
1495
1496IF (this%arraysize > 0) THEN
1497 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1498ELSE
1499 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1500ENDIF
1501IF (shpfileisnull(shphandle)) THEN
1502 ! log here
1503 CALL raise_error()
1504 RETURN
1505ENDIF
1506
1507DO i = 1, this%arraysize
1508 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1509ENDDO
1510
1511CALL shpclose(shphandle)
1512
1513END SUBROUTINE arrayof_georef_coord_array_export
1514#endif
1515
1527FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1528TYPE(georef_coord), INTENT(IN) :: this
1529TYPE(georef_coord_array), INTENT(IN) :: poly
1530LOGICAL :: inside
1531
1532INTEGER :: i
1533
1534inside = .false.
1536IF (.NOT.ALLOCATED(poly%coord)) RETURN
1537! if outside bounding box stop here
1538IF (poly%bbox_updated) THEN
1539 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1540ENDIF
1541
1542IF (ALLOCATED(poly%parts)) THEN
1543 DO i = 1, SIZE(poly%parts)-1
1545 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1546 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1547 ENDDO
1548 IF (SIZE(poly%parts) > 0) THEN ! safety check
1550 poly%coord(poly%parts(i)+1:)%x, &
1551 poly%coord(poly%parts(i)+1:)%y)
1552 ENDIF
1553
1554ELSE
1555 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1556 inside = pointinpoly(this%x, this%y, &
1557 poly%coord(:)%x, poly%coord(:)%y)
1558ENDIF
1559
1560CONTAINS
1561
1562FUNCTION pointinpoly(x, y, px, py)
1563DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1564LOGICAL :: pointinpoly
1565
1566INTEGER :: i, j, starti
1567
1568pointinpoly = .false.
1569
1570IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1571 starti = 2
1572 j = 1
1573ELSE ! unclosed polygon
1574 starti = 1
1575 j = SIZE(px)
1576ENDIF
1577DO i = starti, SIZE(px)
1578 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1579 (py(j) <= y .AND. y < py(i))) THEN
1580 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1581 pointinpoly = .NOT. pointinpoly
1582 ENDIF
1583 ENDIF
1584 j = i
1585ENDDO
1586
1587END FUNCTION pointinpoly
1588
1589END FUNCTION georef_coord_inside
1590
1591
1592
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 |