libsim Versione 7.2.0
|
◆ 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 864 del file georef_coord_class.F90. 865! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
866! authors:
867! Davide Cesari <dcesari@arpa.emr.it>
868! Paolo Patruno <ppatruno@arpa.emr.it>
869
870! This program is free software; you can redistribute it and/or
871! modify it under the terms of the GNU General Public License as
872! published by the Free Software Foundation; either version 2 of
873! the License, or (at your option) any later version.
874
875! This program is distributed in the hope that it will be useful,
876! but WITHOUT ANY WARRANTY; without even the implied warranty of
877! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
878! GNU General Public License for more details.
879
880! You should have received a copy of the GNU General Public License
881! along with this program. If not, see <http://www.gnu.org/licenses/>.
882#include "config.h"
883
901USE geo_proj_class
902#ifdef HAVE_SHAPELIB
903USE shapelib
904#endif
905IMPLICIT NONE
906
912 PRIVATE
913 DOUBLE PRECISION :: x=dmiss, y=dmiss
915
918
924 PRIVATE
925 INTEGER,ALLOCATABLE :: parts(:)
926 TYPE(georef_coord),ALLOCATABLE :: coord(:)
927 INTEGER :: topo=imiss
928 TYPE(geo_proj) :: proj
929 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
930 LOGICAL :: bbox_updated=.false.
932
933INTEGER,PARAMETER :: georef_coord_array_point = 1
934INTEGER,PARAMETER :: georef_coord_array_arc = 3
935INTEGER,PARAMETER :: georef_coord_array_polygon = 5
936INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
937
938
943 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
944END INTERFACE
945
948 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
949END INTERFACE
950
953 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
954END INTERFACE
955
956INTERFACE compute_bbox
957 MODULE PROCEDURE georef_coord_array_compute_bbox
958END INTERFACE
959
961INTERFACE OPERATOR (==)
962 MODULE PROCEDURE georef_coord_eq
963END INTERFACE
964
966INTERFACE OPERATOR (/=)
967 MODULE PROCEDURE georef_coord_ne
968END INTERFACE
969
972INTERFACE OPERATOR (>=)
973 MODULE PROCEDURE georef_coord_ge
974END INTERFACE
975
978INTERFACE OPERATOR (<=)
979 MODULE PROCEDURE georef_coord_le
980END INTERFACE
981
982#ifdef HAVE_SHAPELIB
983
985INTERFACE import
986 MODULE PROCEDURE arrayof_georef_coord_array_import
987END INTERFACE
988
992 MODULE PROCEDURE arrayof_georef_coord_array_export
993END INTERFACE
994#endif
995
999 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1000END INTERFACE
1001
1005 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1006END INTERFACE
1007
1010 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1011END INTERFACE
1012
1015 MODULE PROCEDURE georef_coord_dist
1016END INTERFACE
1017
1018#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1019#define ARRAYOF_TYPE arrayof_georef_coord_array
1020!define ARRAYOF_ORIGEQ 0
1021#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1022#include "arrayof_pre.F90"
1023! from arrayof
1025
1026PRIVATE
1028 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1029 georef_coord_array_polygon, georef_coord_array_multipoint, &
1031#ifdef HAVE_SHAPELIB
1033#endif
1035 georef_coord_new, georef_coord_array_new
1036
1037CONTAINS
1038
1039#include "arrayof_post.F90"
1040
1041! ===================
1042! == georef_coord ==
1043! ===================
1047FUNCTION georef_coord_new(x, y) RESULT(this)
1048DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1049DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1050TYPE(georef_coord) :: this
1051
1054
1055END FUNCTION georef_coord_new
1056
1057
1058SUBROUTINE georef_coord_delete(this)
1059TYPE(georef_coord),INTENT(inout) :: this
1060
1061this%x = dmiss
1062this%y = dmiss
1063
1064END SUBROUTINE georef_coord_delete
1065
1066
1067ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1068TYPE(georef_coord),INTENT(in) :: this
1069LOGICAL :: res
1070
1071res = .NOT. this == georef_coord_miss
1072
1073END FUNCTION georef_coord_c_e
1074
1075
1082ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1083TYPE(georef_coord),INTENT(in) :: this
1084DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1085DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1086
1087IF (PRESENT(x)) x = this%x
1088IF (PRESENT(y)) y = this%y
1089
1090END SUBROUTINE georef_coord_getval
1091
1092
1101ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1102TYPE(georef_coord),INTENT(in) :: this
1103TYPE(geo_proj),INTENT(in) :: proj
1104DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1105DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1106DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1107DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1108
1109DOUBLE PRECISION :: llon, llat
1110
1111IF (PRESENT(x)) x = this%x
1112IF (PRESENT(y)) y = this%y
1113IF (PRESENT(lon) .OR. present(lat)) THEN
1115 IF (PRESENT(lon)) lon = llon
1116 IF (PRESENT(lat)) lat = llat
1117ENDIF
1118
1119END SUBROUTINE georef_coord_proj_getval
1120
1121
1122! document and improve
1123ELEMENTAL FUNCTION getlat(this)
1124TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1125DOUBLE PRECISION :: getlat ! latitudine geografica
1126
1127getlat = this%y ! change!!!
1128
1129END FUNCTION getlat
1130
1131! document and improve
1132ELEMENTAL FUNCTION getlon(this)
1133TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1134DOUBLE PRECISION :: getlon ! longitudine geografica
1135
1136getlon = this%x ! change!!!
1137
1138END FUNCTION getlon
1139
1140
1141ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1142TYPE(georef_coord),INTENT(IN) :: this, that
1143LOGICAL :: res
1144
1145res = (this%x == that%x .AND. this%y == that%y)
1146
1147END FUNCTION georef_coord_eq
1148
1149
1150ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1151TYPE(georef_coord),INTENT(IN) :: this, that
1152LOGICAL :: res
1153
1154res = (this%x >= that%x .AND. this%y >= that%y)
1155
1156END FUNCTION georef_coord_ge
1157
1158
1159ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1160TYPE(georef_coord),INTENT(IN) :: this, that
1161LOGICAL :: res
1162
1163res = (this%x <= that%x .AND. this%y <= that%y)
1164
1165END FUNCTION georef_coord_le
1166
1167
1168ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1169TYPE(georef_coord),INTENT(IN) :: this, that
1170LOGICAL :: res
1171
1172res = .NOT.(this == that)
1173
1174END FUNCTION georef_coord_ne
1175
1176
1182SUBROUTINE georef_coord_read_unit(this, unit)
1183TYPE(georef_coord),INTENT(out) :: this
1184INTEGER, INTENT(in) :: unit
1185
1186CALL georef_coord_vect_read_unit((/this/), unit)
1187
1188END SUBROUTINE georef_coord_read_unit
1189
1190
1196SUBROUTINE georef_coord_vect_read_unit(this, unit)
1197TYPE(georef_coord) :: this(:)
1198INTEGER, INTENT(in) :: unit
1199
1200CHARACTER(len=40) :: form
1201INTEGER :: i
1202
1203INQUIRE(unit, form=form)
1204IF (form == 'FORMATTED') THEN
1205 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1206!TODO bug gfortran compiler !
1207!missing values are unredeable when formatted
1208ELSE
1209 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1210ENDIF
1211
1212END SUBROUTINE georef_coord_vect_read_unit
1213
1214
1219SUBROUTINE georef_coord_write_unit(this, unit)
1220TYPE(georef_coord),INTENT(in) :: this
1221INTEGER, INTENT(in) :: unit
1222
1223CALL georef_coord_vect_write_unit((/this/), unit)
1224
1225END SUBROUTINE georef_coord_write_unit
1226
1227
1232SUBROUTINE georef_coord_vect_write_unit(this, unit)
1233TYPE(georef_coord),INTENT(in) :: this(:)
1234INTEGER, INTENT(in) :: unit
1235
1236CHARACTER(len=40) :: form
1237INTEGER :: i
1238
1239INQUIRE(unit, form=form)
1240IF (form == 'FORMATTED') THEN
1241 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1242!TODO bug gfortran compiler !
1243!missing values are unredeable when formatted
1244ELSE
1245 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1246ENDIF
1247
1248END SUBROUTINE georef_coord_vect_write_unit
1249
1250
1253FUNCTION georef_coord_dist(this, that) RESULT(dist)
1255TYPE(georef_coord), INTENT (IN) :: this
1256TYPE(georef_coord), INTENT (IN) :: that
1257DOUBLE PRECISION :: dist
1258
1259DOUBLE PRECISION :: x,y
1260! Distanza approssimata, valida per piccoli angoli
1261
1262x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1263y = (this%y-that%y)
1264dist = sqrt(x**2 + y**2)*degrad*rearth
1265
1266END FUNCTION georef_coord_dist
1267
1268
1274FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1275TYPE(georef_coord),INTENT(IN) :: this
1276TYPE(georef_coord),INTENT(IN) :: coordmin
1277TYPE(georef_coord),INTENT(IN) :: coordmax
1278LOGICAL :: res
1279
1280res = (this >= coordmin .AND. this <= coordmax)
1281
1282END FUNCTION georef_coord_inside_rectang
1283
1284
1285! ========================
1286! == georef_coord_array ==
1287! ========================
1293FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1294DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1295DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1296INTEGER,INTENT(in),OPTIONAL :: topo
1297TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1298TYPE(georef_coord_array) :: this
1299
1300INTEGER :: lsize
1301
1302IF (PRESENT(x) .AND. PRESENT(y)) THEN
1303 lsize = min(SIZE(x), SIZE(y))
1304 ALLOCATE(this%coord(lsize))
1305 this%coord(1:lsize)%x = x(1:lsize)
1306 this%coord(1:lsize)%y = y(1:lsize)
1307ENDIF
1308this%topo = optio_l(topo)
1310
1311END FUNCTION georef_coord_array_new
1312
1313
1314SUBROUTINE georef_coord_array_delete(this)
1315TYPE(georef_coord_array),INTENT(inout) :: this
1316
1317TYPE(georef_coord_array) :: lobj
1318
1319this = lobj
1320
1321END SUBROUTINE georef_coord_array_delete
1322
1323
1324ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1325TYPE(georef_coord_array),INTENT(in) :: this
1326LOGICAL :: res
1327
1328res = ALLOCATED(this%coord)
1329
1330END FUNCTION georef_coord_array_c_e
1331
1332
1337SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1338TYPE(georef_coord_array),INTENT(in) :: this
1339DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1340DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1341! allocatable per vedere di nascosto l'effetto che fa
1342INTEGER,OPTIONAL,INTENT(out) :: topo
1343TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1344
1345
1346IF (PRESENT(x)) THEN
1347 IF (ALLOCATED(this%coord)) THEN
1348 x = this%coord%x
1349 ENDIF
1350ENDIF
1351IF (PRESENT(y)) THEN
1352 IF (ALLOCATED(this%coord)) THEN
1353 y = this%coord%y
1354 ENDIF
1355ENDIF
1356IF (PRESENT(topo)) topo = this%topo
1358
1359END SUBROUTINE georef_coord_array_getval
1360
1361
1367SUBROUTINE georef_coord_array_compute_bbox(this)
1368TYPE(georef_coord_array),INTENT(inout) :: this
1369
1370IF (ALLOCATED(this%coord)) THEN
1371 this%bbox(1)%x = minval(this%coord(:)%x)
1372 this%bbox(1)%y = minval(this%coord(:)%y)
1373 this%bbox(2)%x = maxval(this%coord(:)%x)
1374 this%bbox(2)%y = maxval(this%coord(:)%y)
1375 this%bbox_updated = .true.
1376ENDIF
1377
1378END SUBROUTINE georef_coord_array_compute_bbox
1379
1380#ifdef HAVE_SHAPELIB
1381! internal method for importing a single shape
1382SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1383TYPE(georef_coord_array),INTENT(OUT) :: this
1384TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1385INTEGER,INTENT(IN) :: nshp
1386
1387TYPE(shpobject) :: shpobj
1388
1389! read shape object
1390shpobj = shpreadobject(shphandle, nshp)
1391IF (.NOT.shpisnull(shpobj)) THEN
1392! import it in georef_coord object
1393 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1394 topo=shpobj%nshptype)
1395 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1396 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1397 ELSE IF (ALLOCATED(this%parts)) THEN
1398 DEALLOCATE(this%parts)
1399 ENDIF
1400 CALL shpdestroyobject(shpobj)
1401 CALL compute_bbox(this)
1402ENDIF
1403
1404
1405END SUBROUTINE georef_coord_array_import
1406
1407
1408! internal method for exporting a single shape
1409SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1410TYPE(georef_coord_array),INTENT(in) :: this
1411TYPE(shpfileobject),INTENT(inout) :: shphandle
1412INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1413
1414INTEGER :: i
1415TYPE(shpobject) :: shpobj
1416
1417IF (ALLOCATED(this%coord)) THEN
1418 IF (ALLOCATED(this%parts)) THEN
1419 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1420 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1421 ELSE
1422 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1423 this%coord(:)%x, this%coord(:)%y)
1424 ENDIF
1425ELSE
1426 RETURN
1427ENDIF
1428
1429IF (.NOT.shpisnull(shpobj)) THEN
1430 i = shpwriteobject(shphandle, nshp, shpobj)
1431 CALL shpdestroyobject(shpobj)
1432ENDIF
1433
1434END SUBROUTINE georef_coord_array_export
1435
1436
1447SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1448TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1449CHARACTER(len=*),INTENT(in) :: shpfile
1450
1451REAL(kind=fp_d) :: minb(4), maxb(4)
1452INTEGER :: i, ns, shptype, dbfnf, dbfnr
1453TYPE(shpfileobject) :: shphandle
1454
1455shphandle = shpopen(trim(shpfile), 'rb')
1456IF (shpfileisnull(shphandle)) THEN
1457 ! log here
1458 CALL raise_error()
1459 RETURN
1460ENDIF
1461
1462! get info about file
1463CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1464IF (ns > 0) THEN ! allocate and read the object
1466 DO i = 1, ns
1467 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1468 ENDDO
1469ENDIF
1470
1471CALL shpclose(shphandle)
1472! pack object to save memory
1474
1475END SUBROUTINE arrayof_georef_coord_array_import
1476
1477
1483SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1484TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1485CHARACTER(len=*),INTENT(in) :: shpfile
1486
1487INTEGER :: i
1488TYPE(shpfileobject) :: shphandle
1489
1490IF (this%arraysize > 0) THEN
1491 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1492ELSE
1493 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1494ENDIF
1495IF (shpfileisnull(shphandle)) THEN
1496 ! log here
1497 CALL raise_error()
1498 RETURN
1499ENDIF
1500
1501DO i = 1, this%arraysize
1502 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1503ENDDO
1504
1505CALL shpclose(shphandle)
1506
1507END SUBROUTINE arrayof_georef_coord_array_export
1508#endif
1509
1521FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1522TYPE(georef_coord), INTENT(IN) :: this
1523TYPE(georef_coord_array), INTENT(IN) :: poly
1524LOGICAL :: inside
1525
1526INTEGER :: i
1527
1528inside = .false.
1530IF (.NOT.ALLOCATED(poly%coord)) RETURN
1531! if outside bounding box stop here
1532IF (poly%bbox_updated) THEN
1533 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1534ENDIF
1535
1536IF (ALLOCATED(poly%parts)) THEN
1537 DO i = 1, SIZE(poly%parts)-1
1539 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1540 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1541 ENDDO
1542 IF (SIZE(poly%parts) > 0) THEN ! safety check
1544 poly%coord(poly%parts(i)+1:)%x, &
1545 poly%coord(poly%parts(i)+1:)%y)
1546 ENDIF
1547
1548ELSE
1549 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1550 inside = pointinpoly(this%x, this%y, &
1551 poly%coord(:)%x, poly%coord(:)%y)
1552ENDIF
1553
1554CONTAINS
1555
1556FUNCTION pointinpoly(x, y, px, py)
1557DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1558LOGICAL :: pointinpoly
1559
1560INTEGER :: i, j, starti
1561
1562pointinpoly = .false.
1563
1564IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1565 starti = 2
1566 j = 1
1567ELSE ! unclosed polygon
1568 starti = 1
1569 j = SIZE(px)
1570ENDIF
1571DO i = starti, SIZE(px)
1572 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1573 (py(j) <= y .AND. y < py(i))) THEN
1574 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1575 pointinpoly = .NOT. pointinpoly
1576 ENDIF
1577 ENDIF
1578 j = i
1579ENDDO
1580
1581END FUNCTION pointinpoly
1582
1583END FUNCTION georef_coord_inside
1584
1585
1586
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 |