libsim Versione 7.2.4
|
◆ georef_coord_vect_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 900 del file georef_coord_class.F90. 901! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
902! authors:
903! Davide Cesari <dcesari@arpa.emr.it>
904! Paolo Patruno <ppatruno@arpa.emr.it>
905
906! This program is free software; you can redistribute it and/or
907! modify it under the terms of the GNU General Public License as
908! published by the Free Software Foundation; either version 2 of
909! the License, or (at your option) any later version.
910
911! This program is distributed in the hope that it will be useful,
912! but WITHOUT ANY WARRANTY; without even the implied warranty of
913! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
914! GNU General Public License for more details.
915
916! You should have received a copy of the GNU General Public License
917! along with this program. If not, see <http://www.gnu.org/licenses/>.
918#include "config.h"
919
937USE geo_proj_class
938#ifdef HAVE_SHAPELIB
939USE shapelib
940#endif
941IMPLICIT NONE
942
948 PRIVATE
949 DOUBLE PRECISION :: x=dmiss, y=dmiss
951
954
960 PRIVATE
961 INTEGER,ALLOCATABLE :: parts(:)
962 TYPE(georef_coord),ALLOCATABLE :: coord(:)
963 INTEGER :: topo=imiss
964 TYPE(geo_proj) :: proj
965 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
966 LOGICAL :: bbox_updated=.false.
968
969INTEGER,PARAMETER :: georef_coord_array_point = 1
970INTEGER,PARAMETER :: georef_coord_array_arc = 3
971INTEGER,PARAMETER :: georef_coord_array_polygon = 5
972INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
973
974
979 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
980END INTERFACE
981
984 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
985END INTERFACE
986
989 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
990END INTERFACE
991
992INTERFACE compute_bbox
993 MODULE PROCEDURE georef_coord_array_compute_bbox
994END INTERFACE
995
997INTERFACE OPERATOR (==)
998 MODULE PROCEDURE georef_coord_eq
999END INTERFACE
1000
1002INTERFACE OPERATOR (/=)
1003 MODULE PROCEDURE georef_coord_ne
1004END INTERFACE
1005
1008INTERFACE OPERATOR (>=)
1009 MODULE PROCEDURE georef_coord_ge
1010END INTERFACE
1011
1014INTERFACE OPERATOR (<=)
1015 MODULE PROCEDURE georef_coord_le
1016END INTERFACE
1017
1018#ifdef HAVE_SHAPELIB
1019
1021INTERFACE import
1022 MODULE PROCEDURE arrayof_georef_coord_array_import
1023END INTERFACE
1024
1028 MODULE PROCEDURE arrayof_georef_coord_array_export
1029END INTERFACE
1030#endif
1031
1035 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1036END INTERFACE
1037
1041 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1042END INTERFACE
1043
1046 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1047END INTERFACE
1048
1051 MODULE PROCEDURE georef_coord_dist
1052END INTERFACE
1053
1054#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1055#define ARRAYOF_TYPE arrayof_georef_coord_array
1056!define ARRAYOF_ORIGEQ 0
1057#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1058#include "arrayof_pre.F90"
1059! from arrayof
1061
1062PRIVATE
1064 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1065 georef_coord_array_polygon, georef_coord_array_multipoint, &
1067#ifdef HAVE_SHAPELIB
1069#endif
1071 georef_coord_new, georef_coord_array_new
1072
1073CONTAINS
1074
1075#include "arrayof_post.F90"
1076
1077! ===================
1078! == georef_coord ==
1079! ===================
1083FUNCTION georef_coord_new(x, y) RESULT(this)
1084DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1085DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1086TYPE(georef_coord) :: this
1087
1090
1091END FUNCTION georef_coord_new
1092
1093
1094SUBROUTINE georef_coord_delete(this)
1095TYPE(georef_coord),INTENT(inout) :: this
1096
1097this%x = dmiss
1098this%y = dmiss
1099
1100END SUBROUTINE georef_coord_delete
1101
1102
1103ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1104TYPE(georef_coord),INTENT(in) :: this
1105LOGICAL :: res
1106
1107res = .NOT. this == georef_coord_miss
1108
1109END FUNCTION georef_coord_c_e
1110
1111
1118ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1119TYPE(georef_coord),INTENT(in) :: this
1120DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1121DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1122
1123IF (PRESENT(x)) x = this%x
1124IF (PRESENT(y)) y = this%y
1125
1126END SUBROUTINE georef_coord_getval
1127
1128
1137ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1138TYPE(georef_coord),INTENT(in) :: this
1139TYPE(geo_proj),INTENT(in) :: proj
1140DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1141DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1142DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1143DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1144
1145DOUBLE PRECISION :: llon, llat
1146
1147IF (PRESENT(x)) x = this%x
1148IF (PRESENT(y)) y = this%y
1149IF (PRESENT(lon) .OR. present(lat)) THEN
1151 IF (PRESENT(lon)) lon = llon
1152 IF (PRESENT(lat)) lat = llat
1153ENDIF
1154
1155END SUBROUTINE georef_coord_proj_getval
1156
1157
1158! document and improve
1159ELEMENTAL FUNCTION getlat(this)
1160TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1161DOUBLE PRECISION :: getlat ! latitudine geografica
1162
1163getlat = this%y ! change!!!
1164
1165END FUNCTION getlat
1166
1167! document and improve
1168ELEMENTAL FUNCTION getlon(this)
1169TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1170DOUBLE PRECISION :: getlon ! longitudine geografica
1171
1172getlon = this%x ! change!!!
1173
1174END FUNCTION getlon
1175
1176
1177ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1178TYPE(georef_coord),INTENT(IN) :: this, that
1179LOGICAL :: res
1180
1181res = (this%x == that%x .AND. this%y == that%y)
1182
1183END FUNCTION georef_coord_eq
1184
1185
1186ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1187TYPE(georef_coord),INTENT(IN) :: this, that
1188LOGICAL :: res
1189
1190res = (this%x >= that%x .AND. this%y >= that%y)
1191
1192END FUNCTION georef_coord_ge
1193
1194
1195ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1196TYPE(georef_coord),INTENT(IN) :: this, that
1197LOGICAL :: res
1198
1199res = (this%x <= that%x .AND. this%y <= that%y)
1200
1201END FUNCTION georef_coord_le
1202
1203
1204ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1205TYPE(georef_coord),INTENT(IN) :: this, that
1206LOGICAL :: res
1207
1208res = .NOT.(this == that)
1209
1210END FUNCTION georef_coord_ne
1211
1212
1218SUBROUTINE georef_coord_read_unit(this, unit)
1219TYPE(georef_coord),INTENT(out) :: this
1220INTEGER, INTENT(in) :: unit
1221
1222CALL georef_coord_vect_read_unit((/this/), unit)
1223
1224END SUBROUTINE georef_coord_read_unit
1225
1226
1232SUBROUTINE georef_coord_vect_read_unit(this, unit)
1233TYPE(georef_coord) :: this(:)
1234INTEGER, INTENT(in) :: unit
1235
1236CHARACTER(len=40) :: form
1237INTEGER :: i
1238
1239INQUIRE(unit, form=form)
1240IF (form == 'FORMATTED') THEN
1241 read(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 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1246ENDIF
1247
1248END SUBROUTINE georef_coord_vect_read_unit
1249
1250
1255SUBROUTINE georef_coord_write_unit(this, unit)
1256TYPE(georef_coord),INTENT(in) :: this
1257INTEGER, INTENT(in) :: unit
1258
1259CALL georef_coord_vect_write_unit((/this/), unit)
1260
1261END SUBROUTINE georef_coord_write_unit
1262
1263
1268SUBROUTINE georef_coord_vect_write_unit(this, unit)
1269TYPE(georef_coord),INTENT(in) :: this(:)
1270INTEGER, INTENT(in) :: unit
1271
1272CHARACTER(len=40) :: form
1273INTEGER :: i
1274
1275INQUIRE(unit, form=form)
1276IF (form == 'FORMATTED') THEN
1277 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1278!TODO bug gfortran compiler !
1279!missing values are unredeable when formatted
1280ELSE
1281 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1282ENDIF
1283
1284END SUBROUTINE georef_coord_vect_write_unit
1285
1286
1289FUNCTION georef_coord_dist(this, that) RESULT(dist)
1291TYPE(georef_coord), INTENT (IN) :: this
1292TYPE(georef_coord), INTENT (IN) :: that
1293DOUBLE PRECISION :: dist
1294
1295DOUBLE PRECISION :: x,y
1296! Distanza approssimata, valida per piccoli angoli
1297
1298x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1299y = (this%y-that%y)
1300dist = sqrt(x**2 + y**2)*degrad*rearth
1301
1302END FUNCTION georef_coord_dist
1303
1304
1310FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1311TYPE(georef_coord),INTENT(IN) :: this
1312TYPE(georef_coord),INTENT(IN) :: coordmin
1313TYPE(georef_coord),INTENT(IN) :: coordmax
1314LOGICAL :: res
1315
1316res = (this >= coordmin .AND. this <= coordmax)
1317
1318END FUNCTION georef_coord_inside_rectang
1319
1320
1321! ========================
1322! == georef_coord_array ==
1323! ========================
1329FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1330DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1331DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1332INTEGER,INTENT(in),OPTIONAL :: topo
1333TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1334TYPE(georef_coord_array) :: this
1335
1336INTEGER :: lsize
1337
1338IF (PRESENT(x) .AND. PRESENT(y)) THEN
1339 lsize = min(SIZE(x), SIZE(y))
1340 ALLOCATE(this%coord(lsize))
1341 this%coord(1:lsize)%x = x(1:lsize)
1342 this%coord(1:lsize)%y = y(1:lsize)
1343ENDIF
1344this%topo = optio_l(topo)
1346
1347END FUNCTION georef_coord_array_new
1348
1349
1350SUBROUTINE georef_coord_array_delete(this)
1351TYPE(georef_coord_array),INTENT(inout) :: this
1352
1353TYPE(georef_coord_array) :: lobj
1354
1355this = lobj
1356
1357END SUBROUTINE georef_coord_array_delete
1358
1359
1360ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1361TYPE(georef_coord_array),INTENT(in) :: this
1362LOGICAL :: res
1363
1364res = ALLOCATED(this%coord)
1365
1366END FUNCTION georef_coord_array_c_e
1367
1368
1373SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1374TYPE(georef_coord_array),INTENT(in) :: this
1375DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1376DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1377! allocatable per vedere di nascosto l'effetto che fa
1378INTEGER,OPTIONAL,INTENT(out) :: topo
1379TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1380
1381
1382IF (PRESENT(x)) THEN
1383 IF (ALLOCATED(this%coord)) THEN
1384 x = this%coord%x
1385 ENDIF
1386ENDIF
1387IF (PRESENT(y)) THEN
1388 IF (ALLOCATED(this%coord)) THEN
1389 y = this%coord%y
1390 ENDIF
1391ENDIF
1392IF (PRESENT(topo)) topo = this%topo
1394
1395END SUBROUTINE georef_coord_array_getval
1396
1397
1403SUBROUTINE georef_coord_array_compute_bbox(this)
1404TYPE(georef_coord_array),INTENT(inout) :: this
1405
1406IF (ALLOCATED(this%coord)) THEN
1407 this%bbox(1)%x = minval(this%coord(:)%x)
1408 this%bbox(1)%y = minval(this%coord(:)%y)
1409 this%bbox(2)%x = maxval(this%coord(:)%x)
1410 this%bbox(2)%y = maxval(this%coord(:)%y)
1411 this%bbox_updated = .true.
1412ENDIF
1413
1414END SUBROUTINE georef_coord_array_compute_bbox
1415
1416#ifdef HAVE_SHAPELIB
1417! internal method for importing a single shape
1418SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1419TYPE(georef_coord_array),INTENT(OUT) :: this
1420TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1421INTEGER,INTENT(IN) :: nshp
1422
1423TYPE(shpobject) :: shpobj
1424
1425! read shape object
1426shpobj = shpreadobject(shphandle, nshp)
1427IF (.NOT.shpisnull(shpobj)) THEN
1428! import it in georef_coord object
1429 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1430 topo=shpobj%nshptype)
1431 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1432 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1433 ELSE IF (ALLOCATED(this%parts)) THEN
1434 DEALLOCATE(this%parts)
1435 ENDIF
1436 CALL shpdestroyobject(shpobj)
1437 CALL compute_bbox(this)
1438ENDIF
1439
1440
1441END SUBROUTINE georef_coord_array_import
1442
1443
1444! internal method for exporting a single shape
1445SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1446TYPE(georef_coord_array),INTENT(in) :: this
1447TYPE(shpfileobject),INTENT(inout) :: shphandle
1448INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1449
1450INTEGER :: i
1451TYPE(shpobject) :: shpobj
1452
1453IF (ALLOCATED(this%coord)) THEN
1454 IF (ALLOCATED(this%parts)) THEN
1455 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1456 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1457 ELSE
1458 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1459 this%coord(:)%x, this%coord(:)%y)
1460 ENDIF
1461ELSE
1462 RETURN
1463ENDIF
1464
1465IF (.NOT.shpisnull(shpobj)) THEN
1466 i = shpwriteobject(shphandle, nshp, shpobj)
1467 CALL shpdestroyobject(shpobj)
1468ENDIF
1469
1470END SUBROUTINE georef_coord_array_export
1471
1472
1483SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1484TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1485CHARACTER(len=*),INTENT(in) :: shpfile
1486
1487REAL(kind=fp_d) :: minb(4), maxb(4)
1488INTEGER :: i, ns, shptype, dbfnf, dbfnr
1489TYPE(shpfileobject) :: shphandle
1490
1491shphandle = shpopen(trim(shpfile), 'rb')
1492IF (shpfileisnull(shphandle)) THEN
1493 ! log here
1494 CALL raise_error()
1495 RETURN
1496ENDIF
1497
1498! get info about file
1499CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1500IF (ns > 0) THEN ! allocate and read the object
1502 DO i = 1, ns
1503 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1504 ENDDO
1505ENDIF
1506
1507CALL shpclose(shphandle)
1508! pack object to save memory
1510
1511END SUBROUTINE arrayof_georef_coord_array_import
1512
1513
1519SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1520TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1521CHARACTER(len=*),INTENT(in) :: shpfile
1522
1523INTEGER :: i
1524TYPE(shpfileobject) :: shphandle
1525
1526IF (this%arraysize > 0) THEN
1527 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1528ELSE
1529 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1530ENDIF
1531IF (shpfileisnull(shphandle)) THEN
1532 ! log here
1533 CALL raise_error()
1534 RETURN
1535ENDIF
1536
1537DO i = 1, this%arraysize
1538 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1539ENDDO
1540
1541CALL shpclose(shphandle)
1542
1543END SUBROUTINE arrayof_georef_coord_array_export
1544#endif
1545
1557FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1558TYPE(georef_coord), INTENT(IN) :: this
1559TYPE(georef_coord_array), INTENT(IN) :: poly
1560LOGICAL :: inside
1561
1562INTEGER :: i
1563
1564inside = .false.
1566IF (.NOT.ALLOCATED(poly%coord)) RETURN
1567! if outside bounding box stop here
1568IF (poly%bbox_updated) THEN
1569 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1570ENDIF
1571
1572IF (ALLOCATED(poly%parts)) THEN
1573 DO i = 1, SIZE(poly%parts)-1
1575 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1576 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1577 ENDDO
1578 IF (SIZE(poly%parts) > 0) THEN ! safety check
1580 poly%coord(poly%parts(i)+1:)%x, &
1581 poly%coord(poly%parts(i)+1:)%y)
1582 ENDIF
1583
1584ELSE
1585 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1586 inside = pointinpoly(this%x, this%y, &
1587 poly%coord(:)%x, poly%coord(:)%y)
1588ENDIF
1589
1590CONTAINS
1591
1592FUNCTION pointinpoly(x, y, px, py)
1593DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1594LOGICAL :: pointinpoly
1595
1596INTEGER :: i, j, starti
1597
1598pointinpoly = .false.
1599
1600IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1601 starti = 2
1602 j = 1
1603ELSE ! unclosed polygon
1604 starti = 1
1605 j = SIZE(px)
1606ENDIF
1607DO i = starti, SIZE(px)
1608 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1609 (py(j) <= y .AND. y < py(i))) THEN
1610 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1611 pointinpoly = .NOT. pointinpoly
1612 ENDIF
1613 ENDIF
1614 j = i
1615ENDDO
1616
1617END FUNCTION pointinpoly
1618
1619END FUNCTION georef_coord_inside
1620
1621
1622
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 |