libsim Versione 7.2.1
|
◆ georef_coord_write_unit()
Scrive su un'unità di file il contenuto dell'oggetto this. Il record scritto potrà successivamente essere letto con la ::read_unit. Il metodo controlla se il file è aperto per un I/O formattato o non formattato e fa la cosa giusta.
Definizione alla linea 887 del file georef_coord_class.F90. 888! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
889! authors:
890! Davide Cesari <dcesari@arpa.emr.it>
891! Paolo Patruno <ppatruno@arpa.emr.it>
892
893! This program is free software; you can redistribute it and/or
894! modify it under the terms of the GNU General Public License as
895! published by the Free Software Foundation; either version 2 of
896! the License, or (at your option) any later version.
897
898! This program is distributed in the hope that it will be useful,
899! but WITHOUT ANY WARRANTY; without even the implied warranty of
900! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
901! GNU General Public License for more details.
902
903! You should have received a copy of the GNU General Public License
904! along with this program. If not, see <http://www.gnu.org/licenses/>.
905#include "config.h"
906
924USE geo_proj_class
925#ifdef HAVE_SHAPELIB
926USE shapelib
927#endif
928IMPLICIT NONE
929
935 PRIVATE
936 DOUBLE PRECISION :: x=dmiss, y=dmiss
938
941
947 PRIVATE
948 INTEGER,ALLOCATABLE :: parts(:)
949 TYPE(georef_coord),ALLOCATABLE :: coord(:)
950 INTEGER :: topo=imiss
951 TYPE(geo_proj) :: proj
952 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
953 LOGICAL :: bbox_updated=.false.
955
956INTEGER,PARAMETER :: georef_coord_array_point = 1
957INTEGER,PARAMETER :: georef_coord_array_arc = 3
958INTEGER,PARAMETER :: georef_coord_array_polygon = 5
959INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
960
961
966 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
967END INTERFACE
968
971 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
972END INTERFACE
973
976 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
977END INTERFACE
978
979INTERFACE compute_bbox
980 MODULE PROCEDURE georef_coord_array_compute_bbox
981END INTERFACE
982
984INTERFACE OPERATOR (==)
985 MODULE PROCEDURE georef_coord_eq
986END INTERFACE
987
989INTERFACE OPERATOR (/=)
990 MODULE PROCEDURE georef_coord_ne
991END INTERFACE
992
995INTERFACE OPERATOR (>=)
996 MODULE PROCEDURE georef_coord_ge
997END INTERFACE
998
1001INTERFACE OPERATOR (<=)
1002 MODULE PROCEDURE georef_coord_le
1003END INTERFACE
1004
1005#ifdef HAVE_SHAPELIB
1006
1008INTERFACE import
1009 MODULE PROCEDURE arrayof_georef_coord_array_import
1010END INTERFACE
1011
1015 MODULE PROCEDURE arrayof_georef_coord_array_export
1016END INTERFACE
1017#endif
1018
1022 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1023END INTERFACE
1024
1028 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1029END INTERFACE
1030
1033 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1034END INTERFACE
1035
1038 MODULE PROCEDURE georef_coord_dist
1039END INTERFACE
1040
1041#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1042#define ARRAYOF_TYPE arrayof_georef_coord_array
1043!define ARRAYOF_ORIGEQ 0
1044#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1045#include "arrayof_pre.F90"
1046! from arrayof
1048
1049PRIVATE
1051 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1052 georef_coord_array_polygon, georef_coord_array_multipoint, &
1054#ifdef HAVE_SHAPELIB
1056#endif
1058 georef_coord_new, georef_coord_array_new
1059
1060CONTAINS
1061
1062#include "arrayof_post.F90"
1063
1064! ===================
1065! == georef_coord ==
1066! ===================
1070FUNCTION georef_coord_new(x, y) RESULT(this)
1071DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1072DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1073TYPE(georef_coord) :: this
1074
1077
1078END FUNCTION georef_coord_new
1079
1080
1081SUBROUTINE georef_coord_delete(this)
1082TYPE(georef_coord),INTENT(inout) :: this
1083
1084this%x = dmiss
1085this%y = dmiss
1086
1087END SUBROUTINE georef_coord_delete
1088
1089
1090ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1091TYPE(georef_coord),INTENT(in) :: this
1092LOGICAL :: res
1093
1094res = .NOT. this == georef_coord_miss
1095
1096END FUNCTION georef_coord_c_e
1097
1098
1105ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1106TYPE(georef_coord),INTENT(in) :: this
1107DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1108DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1109
1110IF (PRESENT(x)) x = this%x
1111IF (PRESENT(y)) y = this%y
1112
1113END SUBROUTINE georef_coord_getval
1114
1115
1124ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1125TYPE(georef_coord),INTENT(in) :: this
1126TYPE(geo_proj),INTENT(in) :: proj
1127DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1128DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1129DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1130DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1131
1132DOUBLE PRECISION :: llon, llat
1133
1134IF (PRESENT(x)) x = this%x
1135IF (PRESENT(y)) y = this%y
1136IF (PRESENT(lon) .OR. present(lat)) THEN
1138 IF (PRESENT(lon)) lon = llon
1139 IF (PRESENT(lat)) lat = llat
1140ENDIF
1141
1142END SUBROUTINE georef_coord_proj_getval
1143
1144
1145! document and improve
1146ELEMENTAL FUNCTION getlat(this)
1147TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1148DOUBLE PRECISION :: getlat ! latitudine geografica
1149
1150getlat = this%y ! change!!!
1151
1152END FUNCTION getlat
1153
1154! document and improve
1155ELEMENTAL FUNCTION getlon(this)
1156TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1157DOUBLE PRECISION :: getlon ! longitudine geografica
1158
1159getlon = this%x ! change!!!
1160
1161END FUNCTION getlon
1162
1163
1164ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1165TYPE(georef_coord),INTENT(IN) :: this, that
1166LOGICAL :: res
1167
1168res = (this%x == that%x .AND. this%y == that%y)
1169
1170END FUNCTION georef_coord_eq
1171
1172
1173ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1174TYPE(georef_coord),INTENT(IN) :: this, that
1175LOGICAL :: res
1176
1177res = (this%x >= that%x .AND. this%y >= that%y)
1178
1179END FUNCTION georef_coord_ge
1180
1181
1182ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1183TYPE(georef_coord),INTENT(IN) :: this, that
1184LOGICAL :: res
1185
1186res = (this%x <= that%x .AND. this%y <= that%y)
1187
1188END FUNCTION georef_coord_le
1189
1190
1191ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1192TYPE(georef_coord),INTENT(IN) :: this, that
1193LOGICAL :: res
1194
1195res = .NOT.(this == that)
1196
1197END FUNCTION georef_coord_ne
1198
1199
1205SUBROUTINE georef_coord_read_unit(this, unit)
1206TYPE(georef_coord),INTENT(out) :: this
1207INTEGER, INTENT(in) :: unit
1208
1209CALL georef_coord_vect_read_unit((/this/), unit)
1210
1211END SUBROUTINE georef_coord_read_unit
1212
1213
1219SUBROUTINE georef_coord_vect_read_unit(this, unit)
1220TYPE(georef_coord) :: this(:)
1221INTEGER, INTENT(in) :: unit
1222
1223CHARACTER(len=40) :: form
1224INTEGER :: i
1225
1226INQUIRE(unit, form=form)
1227IF (form == 'FORMATTED') THEN
1228 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1229!TODO bug gfortran compiler !
1230!missing values are unredeable when formatted
1231ELSE
1232 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1233ENDIF
1234
1235END SUBROUTINE georef_coord_vect_read_unit
1236
1237
1242SUBROUTINE georef_coord_write_unit(this, unit)
1243TYPE(georef_coord),INTENT(in) :: this
1244INTEGER, INTENT(in) :: unit
1245
1246CALL georef_coord_vect_write_unit((/this/), unit)
1247
1248END SUBROUTINE georef_coord_write_unit
1249
1250
1255SUBROUTINE georef_coord_vect_write_unit(this, unit)
1256TYPE(georef_coord),INTENT(in) :: this(:)
1257INTEGER, INTENT(in) :: unit
1258
1259CHARACTER(len=40) :: form
1260INTEGER :: i
1261
1262INQUIRE(unit, form=form)
1263IF (form == 'FORMATTED') THEN
1264 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1265!TODO bug gfortran compiler !
1266!missing values are unredeable when formatted
1267ELSE
1268 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1269ENDIF
1270
1271END SUBROUTINE georef_coord_vect_write_unit
1272
1273
1276FUNCTION georef_coord_dist(this, that) RESULT(dist)
1278TYPE(georef_coord), INTENT (IN) :: this
1279TYPE(georef_coord), INTENT (IN) :: that
1280DOUBLE PRECISION :: dist
1281
1282DOUBLE PRECISION :: x,y
1283! Distanza approssimata, valida per piccoli angoli
1284
1285x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1286y = (this%y-that%y)
1287dist = sqrt(x**2 + y**2)*degrad*rearth
1288
1289END FUNCTION georef_coord_dist
1290
1291
1297FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1298TYPE(georef_coord),INTENT(IN) :: this
1299TYPE(georef_coord),INTENT(IN) :: coordmin
1300TYPE(georef_coord),INTENT(IN) :: coordmax
1301LOGICAL :: res
1302
1303res = (this >= coordmin .AND. this <= coordmax)
1304
1305END FUNCTION georef_coord_inside_rectang
1306
1307
1308! ========================
1309! == georef_coord_array ==
1310! ========================
1316FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1317DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1318DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1319INTEGER,INTENT(in),OPTIONAL :: topo
1320TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1321TYPE(georef_coord_array) :: this
1322
1323INTEGER :: lsize
1324
1325IF (PRESENT(x) .AND. PRESENT(y)) THEN
1326 lsize = min(SIZE(x), SIZE(y))
1327 ALLOCATE(this%coord(lsize))
1328 this%coord(1:lsize)%x = x(1:lsize)
1329 this%coord(1:lsize)%y = y(1:lsize)
1330ENDIF
1331this%topo = optio_l(topo)
1333
1334END FUNCTION georef_coord_array_new
1335
1336
1337SUBROUTINE georef_coord_array_delete(this)
1338TYPE(georef_coord_array),INTENT(inout) :: this
1339
1340TYPE(georef_coord_array) :: lobj
1341
1342this = lobj
1343
1344END SUBROUTINE georef_coord_array_delete
1345
1346
1347ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1348TYPE(georef_coord_array),INTENT(in) :: this
1349LOGICAL :: res
1350
1351res = ALLOCATED(this%coord)
1352
1353END FUNCTION georef_coord_array_c_e
1354
1355
1360SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1361TYPE(georef_coord_array),INTENT(in) :: this
1362DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1363DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1364! allocatable per vedere di nascosto l'effetto che fa
1365INTEGER,OPTIONAL,INTENT(out) :: topo
1366TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1367
1368
1369IF (PRESENT(x)) THEN
1370 IF (ALLOCATED(this%coord)) THEN
1371 x = this%coord%x
1372 ENDIF
1373ENDIF
1374IF (PRESENT(y)) THEN
1375 IF (ALLOCATED(this%coord)) THEN
1376 y = this%coord%y
1377 ENDIF
1378ENDIF
1379IF (PRESENT(topo)) topo = this%topo
1381
1382END SUBROUTINE georef_coord_array_getval
1383
1384
1390SUBROUTINE georef_coord_array_compute_bbox(this)
1391TYPE(georef_coord_array),INTENT(inout) :: this
1392
1393IF (ALLOCATED(this%coord)) THEN
1394 this%bbox(1)%x = minval(this%coord(:)%x)
1395 this%bbox(1)%y = minval(this%coord(:)%y)
1396 this%bbox(2)%x = maxval(this%coord(:)%x)
1397 this%bbox(2)%y = maxval(this%coord(:)%y)
1398 this%bbox_updated = .true.
1399ENDIF
1400
1401END SUBROUTINE georef_coord_array_compute_bbox
1402
1403#ifdef HAVE_SHAPELIB
1404! internal method for importing a single shape
1405SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1406TYPE(georef_coord_array),INTENT(OUT) :: this
1407TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1408INTEGER,INTENT(IN) :: nshp
1409
1410TYPE(shpobject) :: shpobj
1411
1412! read shape object
1413shpobj = shpreadobject(shphandle, nshp)
1414IF (.NOT.shpisnull(shpobj)) THEN
1415! import it in georef_coord object
1416 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1417 topo=shpobj%nshptype)
1418 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1419 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1420 ELSE IF (ALLOCATED(this%parts)) THEN
1421 DEALLOCATE(this%parts)
1422 ENDIF
1423 CALL shpdestroyobject(shpobj)
1424 CALL compute_bbox(this)
1425ENDIF
1426
1427
1428END SUBROUTINE georef_coord_array_import
1429
1430
1431! internal method for exporting a single shape
1432SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1433TYPE(georef_coord_array),INTENT(in) :: this
1434TYPE(shpfileobject),INTENT(inout) :: shphandle
1435INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1436
1437INTEGER :: i
1438TYPE(shpobject) :: shpobj
1439
1440IF (ALLOCATED(this%coord)) THEN
1441 IF (ALLOCATED(this%parts)) THEN
1442 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1443 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1444 ELSE
1445 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1446 this%coord(:)%x, this%coord(:)%y)
1447 ENDIF
1448ELSE
1449 RETURN
1450ENDIF
1451
1452IF (.NOT.shpisnull(shpobj)) THEN
1453 i = shpwriteobject(shphandle, nshp, shpobj)
1454 CALL shpdestroyobject(shpobj)
1455ENDIF
1456
1457END SUBROUTINE georef_coord_array_export
1458
1459
1470SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1471TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1472CHARACTER(len=*),INTENT(in) :: shpfile
1473
1474REAL(kind=fp_d) :: minb(4), maxb(4)
1475INTEGER :: i, ns, shptype, dbfnf, dbfnr
1476TYPE(shpfileobject) :: shphandle
1477
1478shphandle = shpopen(trim(shpfile), 'rb')
1479IF (shpfileisnull(shphandle)) THEN
1480 ! log here
1481 CALL raise_error()
1482 RETURN
1483ENDIF
1484
1485! get info about file
1486CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1487IF (ns > 0) THEN ! allocate and read the object
1489 DO i = 1, ns
1490 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1491 ENDDO
1492ENDIF
1493
1494CALL shpclose(shphandle)
1495! pack object to save memory
1497
1498END SUBROUTINE arrayof_georef_coord_array_import
1499
1500
1506SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1507TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1508CHARACTER(len=*),INTENT(in) :: shpfile
1509
1510INTEGER :: i
1511TYPE(shpfileobject) :: shphandle
1512
1513IF (this%arraysize > 0) THEN
1514 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1515ELSE
1516 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1517ENDIF
1518IF (shpfileisnull(shphandle)) THEN
1519 ! log here
1520 CALL raise_error()
1521 RETURN
1522ENDIF
1523
1524DO i = 1, this%arraysize
1525 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1526ENDDO
1527
1528CALL shpclose(shphandle)
1529
1530END SUBROUTINE arrayof_georef_coord_array_export
1531#endif
1532
1544FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1545TYPE(georef_coord), INTENT(IN) :: this
1546TYPE(georef_coord_array), INTENT(IN) :: poly
1547LOGICAL :: inside
1548
1549INTEGER :: i
1550
1551inside = .false.
1553IF (.NOT.ALLOCATED(poly%coord)) RETURN
1554! if outside bounding box stop here
1555IF (poly%bbox_updated) THEN
1556 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1557ENDIF
1558
1559IF (ALLOCATED(poly%parts)) THEN
1560 DO i = 1, SIZE(poly%parts)-1
1562 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1563 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1564 ENDDO
1565 IF (SIZE(poly%parts) > 0) THEN ! safety check
1567 poly%coord(poly%parts(i)+1:)%x, &
1568 poly%coord(poly%parts(i)+1:)%y)
1569 ENDIF
1570
1571ELSE
1572 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1573 inside = pointinpoly(this%x, this%y, &
1574 poly%coord(:)%x, poly%coord(:)%y)
1575ENDIF
1576
1577CONTAINS
1578
1579FUNCTION pointinpoly(x, y, px, py)
1580DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1581LOGICAL :: pointinpoly
1582
1583INTEGER :: i, j, starti
1584
1585pointinpoly = .false.
1586
1587IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1588 starti = 2
1589 j = 1
1590ELSE ! unclosed polygon
1591 starti = 1
1592 j = SIZE(px)
1593ENDIF
1594DO i = starti, SIZE(px)
1595 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1596 (py(j) <= y .AND. y < py(i))) THEN
1597 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1598 pointinpoly = .NOT. pointinpoly
1599 ENDIF
1600 ENDIF
1601 j = i
1602ENDDO
1603
1604END FUNCTION pointinpoly
1605
1606END FUNCTION georef_coord_inside
1607
1608
1609
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 |