libsim Versione 7.1.11
|
◆ 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 893 del file georef_coord_class.F90. 894! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
895! authors:
896! Davide Cesari <dcesari@arpa.emr.it>
897! Paolo Patruno <ppatruno@arpa.emr.it>
898
899! This program is free software; you can redistribute it and/or
900! modify it under the terms of the GNU General Public License as
901! published by the Free Software Foundation; either version 2 of
902! the License, or (at your option) any later version.
903
904! This program is distributed in the hope that it will be useful,
905! but WITHOUT ANY WARRANTY; without even the implied warranty of
906! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
907! GNU General Public License for more details.
908
909! You should have received a copy of the GNU General Public License
910! along with this program. If not, see <http://www.gnu.org/licenses/>.
911#include "config.h"
912
930USE geo_proj_class
931#ifdef HAVE_SHAPELIB
932USE shapelib
933#endif
934IMPLICIT NONE
935
941 PRIVATE
942 DOUBLE PRECISION :: x=dmiss, y=dmiss
944
947
953 PRIVATE
954 INTEGER,ALLOCATABLE :: parts(:)
955 TYPE(georef_coord),ALLOCATABLE :: coord(:)
956 INTEGER :: topo=imiss
957 TYPE(geo_proj) :: proj
958 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
959 LOGICAL :: bbox_updated=.false.
961
962INTEGER,PARAMETER :: georef_coord_array_point = 1
963INTEGER,PARAMETER :: georef_coord_array_arc = 3
964INTEGER,PARAMETER :: georef_coord_array_polygon = 5
965INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
966
967
972 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
973END INTERFACE
974
977 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
978END INTERFACE
979
982 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
983END INTERFACE
984
985INTERFACE compute_bbox
986 MODULE PROCEDURE georef_coord_array_compute_bbox
987END INTERFACE
988
990INTERFACE OPERATOR (==)
991 MODULE PROCEDURE georef_coord_eq
992END INTERFACE
993
995INTERFACE OPERATOR (/=)
996 MODULE PROCEDURE georef_coord_ne
997END INTERFACE
998
1001INTERFACE OPERATOR (>=)
1002 MODULE PROCEDURE georef_coord_ge
1003END INTERFACE
1004
1007INTERFACE OPERATOR (<=)
1008 MODULE PROCEDURE georef_coord_le
1009END INTERFACE
1010
1011#ifdef HAVE_SHAPELIB
1012
1014INTERFACE import
1015 MODULE PROCEDURE arrayof_georef_coord_array_import
1016END INTERFACE
1017
1021 MODULE PROCEDURE arrayof_georef_coord_array_export
1022END INTERFACE
1023#endif
1024
1028 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1029END INTERFACE
1030
1034 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1035END INTERFACE
1036
1039 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1040END INTERFACE
1041
1044 MODULE PROCEDURE georef_coord_dist
1045END INTERFACE
1046
1047#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1048#define ARRAYOF_TYPE arrayof_georef_coord_array
1049!define ARRAYOF_ORIGEQ 0
1050#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1051#include "arrayof_pre.F90"
1052! from arrayof
1054
1055PRIVATE
1057 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1058 georef_coord_array_polygon, georef_coord_array_multipoint, &
1060#ifdef HAVE_SHAPELIB
1062#endif
1064 georef_coord_new, georef_coord_array_new
1065
1066CONTAINS
1067
1068#include "arrayof_post.F90"
1069
1070! ===================
1071! == georef_coord ==
1072! ===================
1076FUNCTION georef_coord_new(x, y) RESULT(this)
1077DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1078DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1079TYPE(georef_coord) :: this
1080
1083
1084END FUNCTION georef_coord_new
1085
1086
1087SUBROUTINE georef_coord_delete(this)
1088TYPE(georef_coord),INTENT(inout) :: this
1089
1090this%x = dmiss
1091this%y = dmiss
1092
1093END SUBROUTINE georef_coord_delete
1094
1095
1096ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1097TYPE(georef_coord),INTENT(in) :: this
1098LOGICAL :: res
1099
1100res = .NOT. this == georef_coord_miss
1101
1102END FUNCTION georef_coord_c_e
1103
1104
1111ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1112TYPE(georef_coord),INTENT(in) :: this
1113DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1114DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1115
1116IF (PRESENT(x)) x = this%x
1117IF (PRESENT(y)) y = this%y
1118
1119END SUBROUTINE georef_coord_getval
1120
1121
1130ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1131TYPE(georef_coord),INTENT(in) :: this
1132TYPE(geo_proj),INTENT(in) :: proj
1133DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1134DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1135DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1136DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1137
1138DOUBLE PRECISION :: llon, llat
1139
1140IF (PRESENT(x)) x = this%x
1141IF (PRESENT(y)) y = this%y
1142IF (PRESENT(lon) .OR. present(lat)) THEN
1144 IF (PRESENT(lon)) lon = llon
1145 IF (PRESENT(lat)) lat = llat
1146ENDIF
1147
1148END SUBROUTINE georef_coord_proj_getval
1149
1150
1151! document and improve
1152ELEMENTAL FUNCTION getlat(this)
1153TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1154DOUBLE PRECISION :: getlat ! latitudine geografica
1155
1156getlat = this%y ! change!!!
1157
1158END FUNCTION getlat
1159
1160! document and improve
1161ELEMENTAL FUNCTION getlon(this)
1162TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1163DOUBLE PRECISION :: getlon ! longitudine geografica
1164
1165getlon = this%x ! change!!!
1166
1167END FUNCTION getlon
1168
1169
1170ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1171TYPE(georef_coord),INTENT(IN) :: this, that
1172LOGICAL :: res
1173
1174res = (this%x == that%x .AND. this%y == that%y)
1175
1176END FUNCTION georef_coord_eq
1177
1178
1179ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1180TYPE(georef_coord),INTENT(IN) :: this, that
1181LOGICAL :: res
1182
1183res = (this%x >= that%x .AND. this%y >= that%y)
1184
1185END FUNCTION georef_coord_ge
1186
1187
1188ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1189TYPE(georef_coord),INTENT(IN) :: this, that
1190LOGICAL :: res
1191
1192res = (this%x <= that%x .AND. this%y <= that%y)
1193
1194END FUNCTION georef_coord_le
1195
1196
1197ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1198TYPE(georef_coord),INTENT(IN) :: this, that
1199LOGICAL :: res
1200
1201res = .NOT.(this == that)
1202
1203END FUNCTION georef_coord_ne
1204
1205
1211SUBROUTINE georef_coord_read_unit(this, unit)
1212TYPE(georef_coord),INTENT(out) :: this
1213INTEGER, INTENT(in) :: unit
1214
1215CALL georef_coord_vect_read_unit((/this/), unit)
1216
1217END SUBROUTINE georef_coord_read_unit
1218
1219
1225SUBROUTINE georef_coord_vect_read_unit(this, unit)
1226TYPE(georef_coord) :: this(:)
1227INTEGER, INTENT(in) :: unit
1228
1229CHARACTER(len=40) :: form
1230INTEGER :: i
1231
1232INQUIRE(unit, form=form)
1233IF (form == 'FORMATTED') THEN
1234 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1235!TODO bug gfortran compiler !
1236!missing values are unredeable when formatted
1237ELSE
1238 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1239ENDIF
1240
1241END SUBROUTINE georef_coord_vect_read_unit
1242
1243
1248SUBROUTINE georef_coord_write_unit(this, unit)
1249TYPE(georef_coord),INTENT(in) :: this
1250INTEGER, INTENT(in) :: unit
1251
1252CALL georef_coord_vect_write_unit((/this/), unit)
1253
1254END SUBROUTINE georef_coord_write_unit
1255
1256
1261SUBROUTINE georef_coord_vect_write_unit(this, unit)
1262TYPE(georef_coord),INTENT(in) :: this(:)
1263INTEGER, INTENT(in) :: unit
1264
1265CHARACTER(len=40) :: form
1266INTEGER :: i
1267
1268INQUIRE(unit, form=form)
1269IF (form == 'FORMATTED') THEN
1270 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1271!TODO bug gfortran compiler !
1272!missing values are unredeable when formatted
1273ELSE
1274 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1275ENDIF
1276
1277END SUBROUTINE georef_coord_vect_write_unit
1278
1279
1282FUNCTION georef_coord_dist(this, that) RESULT(dist)
1284TYPE(georef_coord), INTENT (IN) :: this
1285TYPE(georef_coord), INTENT (IN) :: that
1286DOUBLE PRECISION :: dist
1287
1288DOUBLE PRECISION :: x,y
1289! Distanza approssimata, valida per piccoli angoli
1290
1291x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1292y = (this%y-that%y)
1293dist = sqrt(x**2 + y**2)*degrad*rearth
1294
1295END FUNCTION georef_coord_dist
1296
1297
1303FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1304TYPE(georef_coord),INTENT(IN) :: this
1305TYPE(georef_coord),INTENT(IN) :: coordmin
1306TYPE(georef_coord),INTENT(IN) :: coordmax
1307LOGICAL :: res
1308
1309res = (this >= coordmin .AND. this <= coordmax)
1310
1311END FUNCTION georef_coord_inside_rectang
1312
1313
1314! ========================
1315! == georef_coord_array ==
1316! ========================
1322FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1323DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1324DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1325INTEGER,INTENT(in),OPTIONAL :: topo
1326TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1327TYPE(georef_coord_array) :: this
1328
1329INTEGER :: lsize
1330
1331IF (PRESENT(x) .AND. PRESENT(y)) THEN
1332 lsize = min(SIZE(x), SIZE(y))
1333 ALLOCATE(this%coord(lsize))
1334 this%coord(1:lsize)%x = x(1:lsize)
1335 this%coord(1:lsize)%y = y(1:lsize)
1336ENDIF
1337this%topo = optio_l(topo)
1339
1340END FUNCTION georef_coord_array_new
1341
1342
1343SUBROUTINE georef_coord_array_delete(this)
1344TYPE(georef_coord_array),INTENT(inout) :: this
1345
1346TYPE(georef_coord_array) :: lobj
1347
1348this = lobj
1349
1350END SUBROUTINE georef_coord_array_delete
1351
1352
1353ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1354TYPE(georef_coord_array),INTENT(in) :: this
1355LOGICAL :: res
1356
1357res = ALLOCATED(this%coord)
1358
1359END FUNCTION georef_coord_array_c_e
1360
1361
1366SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1367TYPE(georef_coord_array),INTENT(in) :: this
1368DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1369DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1370! allocatable per vedere di nascosto l'effetto che fa
1371INTEGER,OPTIONAL,INTENT(out) :: topo
1372TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1373
1374
1375IF (PRESENT(x)) THEN
1376 IF (ALLOCATED(this%coord)) THEN
1377 x = this%coord%x
1378 ENDIF
1379ENDIF
1380IF (PRESENT(y)) THEN
1381 IF (ALLOCATED(this%coord)) THEN
1382 y = this%coord%y
1383 ENDIF
1384ENDIF
1385IF (PRESENT(topo)) topo = this%topo
1387
1388END SUBROUTINE georef_coord_array_getval
1389
1390
1396SUBROUTINE georef_coord_array_compute_bbox(this)
1397TYPE(georef_coord_array),INTENT(inout) :: this
1398
1399IF (ALLOCATED(this%coord)) THEN
1400 this%bbox(1)%x = minval(this%coord(:)%x)
1401 this%bbox(1)%y = minval(this%coord(:)%y)
1402 this%bbox(2)%x = maxval(this%coord(:)%x)
1403 this%bbox(2)%y = maxval(this%coord(:)%y)
1404 this%bbox_updated = .true.
1405ENDIF
1406
1407END SUBROUTINE georef_coord_array_compute_bbox
1408
1409#ifdef HAVE_SHAPELIB
1410! internal method for importing a single shape
1411SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1412TYPE(georef_coord_array),INTENT(OUT) :: this
1413TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1414INTEGER,INTENT(IN) :: nshp
1415
1416TYPE(shpobject) :: shpobj
1417
1418! read shape object
1419shpobj = shpreadobject(shphandle, nshp)
1420IF (.NOT.shpisnull(shpobj)) THEN
1421! import it in georef_coord object
1422 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1423 topo=shpobj%nshptype)
1424 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1425 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1426 ELSE IF (ALLOCATED(this%parts)) THEN
1427 DEALLOCATE(this%parts)
1428 ENDIF
1429 CALL shpdestroyobject(shpobj)
1430 CALL compute_bbox(this)
1431ENDIF
1432
1433
1434END SUBROUTINE georef_coord_array_import
1435
1436
1437! internal method for exporting a single shape
1438SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1439TYPE(georef_coord_array),INTENT(in) :: this
1440TYPE(shpfileobject),INTENT(inout) :: shphandle
1441INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1442
1443INTEGER :: i
1444TYPE(shpobject) :: shpobj
1445
1446IF (ALLOCATED(this%coord)) THEN
1447 IF (ALLOCATED(this%parts)) THEN
1448 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1449 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1450 ELSE
1451 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1452 this%coord(:)%x, this%coord(:)%y)
1453 ENDIF
1454ELSE
1455 RETURN
1456ENDIF
1457
1458IF (.NOT.shpisnull(shpobj)) THEN
1459 i = shpwriteobject(shphandle, nshp, shpobj)
1460 CALL shpdestroyobject(shpobj)
1461ENDIF
1462
1463END SUBROUTINE georef_coord_array_export
1464
1465
1476SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1477TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1478CHARACTER(len=*),INTENT(in) :: shpfile
1479
1480REAL(kind=fp_d) :: minb(4), maxb(4)
1481INTEGER :: i, ns, shptype, dbfnf, dbfnr
1482TYPE(shpfileobject) :: shphandle
1483
1484shphandle = shpopen(trim(shpfile), 'rb')
1485IF (shpfileisnull(shphandle)) THEN
1486 ! log here
1487 CALL raise_error()
1488 RETURN
1489ENDIF
1490
1491! get info about file
1492CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1493IF (ns > 0) THEN ! allocate and read the object
1495 DO i = 1, ns
1496 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1497 ENDDO
1498ENDIF
1499
1500CALL shpclose(shphandle)
1501! pack object to save memory
1503
1504END SUBROUTINE arrayof_georef_coord_array_import
1505
1506
1512SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1513TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1514CHARACTER(len=*),INTENT(in) :: shpfile
1515
1516INTEGER :: i
1517TYPE(shpfileobject) :: shphandle
1518
1519IF (this%arraysize > 0) THEN
1520 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1521ELSE
1522 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1523ENDIF
1524IF (shpfileisnull(shphandle)) THEN
1525 ! log here
1526 CALL raise_error()
1527 RETURN
1528ENDIF
1529
1530DO i = 1, this%arraysize
1531 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1532ENDDO
1533
1534CALL shpclose(shphandle)
1535
1536END SUBROUTINE arrayof_georef_coord_array_export
1537#endif
1538
1550FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1551TYPE(georef_coord), INTENT(IN) :: this
1552TYPE(georef_coord_array), INTENT(IN) :: poly
1553LOGICAL :: inside
1554
1555INTEGER :: i
1556
1557inside = .false.
1559IF (.NOT.ALLOCATED(poly%coord)) RETURN
1560! if outside bounding box stop here
1561IF (poly%bbox_updated) THEN
1562 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1563ENDIF
1564
1565IF (ALLOCATED(poly%parts)) THEN
1566 DO i = 1, SIZE(poly%parts)-1
1568 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1569 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1570 ENDDO
1571 IF (SIZE(poly%parts) > 0) THEN ! safety check
1573 poly%coord(poly%parts(i)+1:)%x, &
1574 poly%coord(poly%parts(i)+1:)%y)
1575 ENDIF
1576
1577ELSE
1578 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1579 inside = pointinpoly(this%x, this%y, &
1580 poly%coord(:)%x, poly%coord(:)%y)
1581ENDIF
1582
1583CONTAINS
1584
1585FUNCTION pointinpoly(x, y, px, py)
1586DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1587LOGICAL :: pointinpoly
1588
1589INTEGER :: i, j, starti
1590
1591pointinpoly = .false.
1592
1593IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1594 starti = 2
1595 j = 1
1596ELSE ! unclosed polygon
1597 starti = 1
1598 j = SIZE(px)
1599ENDIF
1600DO i = starti, SIZE(px)
1601 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1602 (py(j) <= y .AND. y < py(i))) THEN
1603 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1604 pointinpoly = .NOT. pointinpoly
1605 ENDIF
1606 ENDIF
1607 j = i
1608ENDDO
1609
1610END FUNCTION pointinpoly
1611
1612END FUNCTION georef_coord_inside
1613
1614
1615
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 |