libsim Versione 7.1.11
|
◆ georef_coord_array_new()
Construct a georef_coord_array object with the optional parameters provided. If coordinates are not provided the object obtained is empty (missing, see c_e function), if coordinate arrays are of different lengths the georef_coord_array is initialised to the shortest length provided.
Definizione alla linea 967 del file georef_coord_class.F90. 968! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
969! authors:
970! Davide Cesari <dcesari@arpa.emr.it>
971! Paolo Patruno <ppatruno@arpa.emr.it>
972
973! This program is free software; you can redistribute it and/or
974! modify it under the terms of the GNU General Public License as
975! published by the Free Software Foundation; either version 2 of
976! the License, or (at your option) any later version.
977
978! This program is distributed in the hope that it will be useful,
979! but WITHOUT ANY WARRANTY; without even the implied warranty of
980! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
981! GNU General Public License for more details.
982
983! You should have received a copy of the GNU General Public License
984! along with this program. If not, see <http://www.gnu.org/licenses/>.
985#include "config.h"
986
1004USE geo_proj_class
1005#ifdef HAVE_SHAPELIB
1006USE shapelib
1007#endif
1008IMPLICIT NONE
1009
1015 PRIVATE
1016 DOUBLE PRECISION :: x=dmiss, y=dmiss
1018
1021
1027 PRIVATE
1028 INTEGER,ALLOCATABLE :: parts(:)
1029 TYPE(georef_coord),ALLOCATABLE :: coord(:)
1030 INTEGER :: topo=imiss
1031 TYPE(geo_proj) :: proj
1032 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
1033 LOGICAL :: bbox_updated=.false.
1035
1036INTEGER,PARAMETER :: georef_coord_array_point = 1
1037INTEGER,PARAMETER :: georef_coord_array_arc = 3
1038INTEGER,PARAMETER :: georef_coord_array_polygon = 5
1039INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
1040
1041
1046 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
1047END INTERFACE
1048
1051 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
1052END INTERFACE
1053
1056 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
1057END INTERFACE
1058
1059INTERFACE compute_bbox
1060 MODULE PROCEDURE georef_coord_array_compute_bbox
1061END INTERFACE
1062
1064INTERFACE OPERATOR (==)
1065 MODULE PROCEDURE georef_coord_eq
1066END INTERFACE
1067
1069INTERFACE OPERATOR (/=)
1070 MODULE PROCEDURE georef_coord_ne
1071END INTERFACE
1072
1075INTERFACE OPERATOR (>=)
1076 MODULE PROCEDURE georef_coord_ge
1077END INTERFACE
1078
1081INTERFACE OPERATOR (<=)
1082 MODULE PROCEDURE georef_coord_le
1083END INTERFACE
1084
1085#ifdef HAVE_SHAPELIB
1086
1088INTERFACE import
1089 MODULE PROCEDURE arrayof_georef_coord_array_import
1090END INTERFACE
1091
1095 MODULE PROCEDURE arrayof_georef_coord_array_export
1096END INTERFACE
1097#endif
1098
1102 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1103END INTERFACE
1104
1108 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1109END INTERFACE
1110
1113 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1114END INTERFACE
1115
1118 MODULE PROCEDURE georef_coord_dist
1119END INTERFACE
1120
1121#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1122#define ARRAYOF_TYPE arrayof_georef_coord_array
1123!define ARRAYOF_ORIGEQ 0
1124#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1125#include "arrayof_pre.F90"
1126! from arrayof
1128
1129PRIVATE
1131 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1132 georef_coord_array_polygon, georef_coord_array_multipoint, &
1134#ifdef HAVE_SHAPELIB
1136#endif
1138 georef_coord_new, georef_coord_array_new
1139
1140CONTAINS
1141
1142#include "arrayof_post.F90"
1143
1144! ===================
1145! == georef_coord ==
1146! ===================
1150FUNCTION georef_coord_new(x, y) RESULT(this)
1151DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1152DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1153TYPE(georef_coord) :: this
1154
1157
1158END FUNCTION georef_coord_new
1159
1160
1161SUBROUTINE georef_coord_delete(this)
1162TYPE(georef_coord),INTENT(inout) :: this
1163
1164this%x = dmiss
1165this%y = dmiss
1166
1167END SUBROUTINE georef_coord_delete
1168
1169
1170ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1171TYPE(georef_coord),INTENT(in) :: this
1172LOGICAL :: res
1173
1174res = .NOT. this == georef_coord_miss
1175
1176END FUNCTION georef_coord_c_e
1177
1178
1185ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1186TYPE(georef_coord),INTENT(in) :: this
1187DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1188DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1189
1190IF (PRESENT(x)) x = this%x
1191IF (PRESENT(y)) y = this%y
1192
1193END SUBROUTINE georef_coord_getval
1194
1195
1204ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1205TYPE(georef_coord),INTENT(in) :: this
1206TYPE(geo_proj),INTENT(in) :: proj
1207DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1208DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1209DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1210DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1211
1212DOUBLE PRECISION :: llon, llat
1213
1214IF (PRESENT(x)) x = this%x
1215IF (PRESENT(y)) y = this%y
1216IF (PRESENT(lon) .OR. present(lat)) THEN
1218 IF (PRESENT(lon)) lon = llon
1219 IF (PRESENT(lat)) lat = llat
1220ENDIF
1221
1222END SUBROUTINE georef_coord_proj_getval
1223
1224
1225! document and improve
1226ELEMENTAL FUNCTION getlat(this)
1227TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1228DOUBLE PRECISION :: getlat ! latitudine geografica
1229
1230getlat = this%y ! change!!!
1231
1232END FUNCTION getlat
1233
1234! document and improve
1235ELEMENTAL FUNCTION getlon(this)
1236TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1237DOUBLE PRECISION :: getlon ! longitudine geografica
1238
1239getlon = this%x ! change!!!
1240
1241END FUNCTION getlon
1242
1243
1244ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1245TYPE(georef_coord),INTENT(IN) :: this, that
1246LOGICAL :: res
1247
1248res = (this%x == that%x .AND. this%y == that%y)
1249
1250END FUNCTION georef_coord_eq
1251
1252
1253ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1254TYPE(georef_coord),INTENT(IN) :: this, that
1255LOGICAL :: res
1256
1257res = (this%x >= that%x .AND. this%y >= that%y)
1258
1259END FUNCTION georef_coord_ge
1260
1261
1262ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1263TYPE(georef_coord),INTENT(IN) :: this, that
1264LOGICAL :: res
1265
1266res = (this%x <= that%x .AND. this%y <= that%y)
1267
1268END FUNCTION georef_coord_le
1269
1270
1271ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1272TYPE(georef_coord),INTENT(IN) :: this, that
1273LOGICAL :: res
1274
1275res = .NOT.(this == that)
1276
1277END FUNCTION georef_coord_ne
1278
1279
1285SUBROUTINE georef_coord_read_unit(this, unit)
1286TYPE(georef_coord),INTENT(out) :: this
1287INTEGER, INTENT(in) :: unit
1288
1289CALL georef_coord_vect_read_unit((/this/), unit)
1290
1291END SUBROUTINE georef_coord_read_unit
1292
1293
1299SUBROUTINE georef_coord_vect_read_unit(this, unit)
1300TYPE(georef_coord) :: this(:)
1301INTEGER, INTENT(in) :: unit
1302
1303CHARACTER(len=40) :: form
1304INTEGER :: i
1305
1306INQUIRE(unit, form=form)
1307IF (form == 'FORMATTED') THEN
1308 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1309!TODO bug gfortran compiler !
1310!missing values are unredeable when formatted
1311ELSE
1312 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1313ENDIF
1314
1315END SUBROUTINE georef_coord_vect_read_unit
1316
1317
1322SUBROUTINE georef_coord_write_unit(this, unit)
1323TYPE(georef_coord),INTENT(in) :: this
1324INTEGER, INTENT(in) :: unit
1325
1326CALL georef_coord_vect_write_unit((/this/), unit)
1327
1328END SUBROUTINE georef_coord_write_unit
1329
1330
1335SUBROUTINE georef_coord_vect_write_unit(this, unit)
1336TYPE(georef_coord),INTENT(in) :: this(:)
1337INTEGER, INTENT(in) :: unit
1338
1339CHARACTER(len=40) :: form
1340INTEGER :: i
1341
1342INQUIRE(unit, form=form)
1343IF (form == 'FORMATTED') THEN
1344 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1345!TODO bug gfortran compiler !
1346!missing values are unredeable when formatted
1347ELSE
1348 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1349ENDIF
1350
1351END SUBROUTINE georef_coord_vect_write_unit
1352
1353
1356FUNCTION georef_coord_dist(this, that) RESULT(dist)
1358TYPE(georef_coord), INTENT (IN) :: this
1359TYPE(georef_coord), INTENT (IN) :: that
1360DOUBLE PRECISION :: dist
1361
1362DOUBLE PRECISION :: x,y
1363! Distanza approssimata, valida per piccoli angoli
1364
1365x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1366y = (this%y-that%y)
1367dist = sqrt(x**2 + y**2)*degrad*rearth
1368
1369END FUNCTION georef_coord_dist
1370
1371
1377FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1378TYPE(georef_coord),INTENT(IN) :: this
1379TYPE(georef_coord),INTENT(IN) :: coordmin
1380TYPE(georef_coord),INTENT(IN) :: coordmax
1381LOGICAL :: res
1382
1383res = (this >= coordmin .AND. this <= coordmax)
1384
1385END FUNCTION georef_coord_inside_rectang
1386
1387
1388! ========================
1389! == georef_coord_array ==
1390! ========================
1396FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1397DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1398DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1399INTEGER,INTENT(in),OPTIONAL :: topo
1400TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1401TYPE(georef_coord_array) :: this
1402
1403INTEGER :: lsize
1404
1405IF (PRESENT(x) .AND. PRESENT(y)) THEN
1406 lsize = min(SIZE(x), SIZE(y))
1407 ALLOCATE(this%coord(lsize))
1408 this%coord(1:lsize)%x = x(1:lsize)
1409 this%coord(1:lsize)%y = y(1:lsize)
1410ENDIF
1411this%topo = optio_l(topo)
1413
1414END FUNCTION georef_coord_array_new
1415
1416
1417SUBROUTINE georef_coord_array_delete(this)
1418TYPE(georef_coord_array),INTENT(inout) :: this
1419
1420TYPE(georef_coord_array) :: lobj
1421
1422this = lobj
1423
1424END SUBROUTINE georef_coord_array_delete
1425
1426
1427ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1428TYPE(georef_coord_array),INTENT(in) :: this
1429LOGICAL :: res
1430
1431res = ALLOCATED(this%coord)
1432
1433END FUNCTION georef_coord_array_c_e
1434
1435
1440SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1441TYPE(georef_coord_array),INTENT(in) :: this
1442DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1443DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1444! allocatable per vedere di nascosto l'effetto che fa
1445INTEGER,OPTIONAL,INTENT(out) :: topo
1446TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1447
1448
1449IF (PRESENT(x)) THEN
1450 IF (ALLOCATED(this%coord)) THEN
1451 x = this%coord%x
1452 ENDIF
1453ENDIF
1454IF (PRESENT(y)) THEN
1455 IF (ALLOCATED(this%coord)) THEN
1456 y = this%coord%y
1457 ENDIF
1458ENDIF
1459IF (PRESENT(topo)) topo = this%topo
1461
1462END SUBROUTINE georef_coord_array_getval
1463
1464
1470SUBROUTINE georef_coord_array_compute_bbox(this)
1471TYPE(georef_coord_array),INTENT(inout) :: this
1472
1473IF (ALLOCATED(this%coord)) THEN
1474 this%bbox(1)%x = minval(this%coord(:)%x)
1475 this%bbox(1)%y = minval(this%coord(:)%y)
1476 this%bbox(2)%x = maxval(this%coord(:)%x)
1477 this%bbox(2)%y = maxval(this%coord(:)%y)
1478 this%bbox_updated = .true.
1479ENDIF
1480
1481END SUBROUTINE georef_coord_array_compute_bbox
1482
1483#ifdef HAVE_SHAPELIB
1484! internal method for importing a single shape
1485SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1486TYPE(georef_coord_array),INTENT(OUT) :: this
1487TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1488INTEGER,INTENT(IN) :: nshp
1489
1490TYPE(shpobject) :: shpobj
1491
1492! read shape object
1493shpobj = shpreadobject(shphandle, nshp)
1494IF (.NOT.shpisnull(shpobj)) THEN
1495! import it in georef_coord object
1496 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1497 topo=shpobj%nshptype)
1498 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1499 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1500 ELSE IF (ALLOCATED(this%parts)) THEN
1501 DEALLOCATE(this%parts)
1502 ENDIF
1503 CALL shpdestroyobject(shpobj)
1504 CALL compute_bbox(this)
1505ENDIF
1506
1507
1508END SUBROUTINE georef_coord_array_import
1509
1510
1511! internal method for exporting a single shape
1512SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1513TYPE(georef_coord_array),INTENT(in) :: this
1514TYPE(shpfileobject),INTENT(inout) :: shphandle
1515INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1516
1517INTEGER :: i
1518TYPE(shpobject) :: shpobj
1519
1520IF (ALLOCATED(this%coord)) THEN
1521 IF (ALLOCATED(this%parts)) THEN
1522 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1523 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1524 ELSE
1525 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1526 this%coord(:)%x, this%coord(:)%y)
1527 ENDIF
1528ELSE
1529 RETURN
1530ENDIF
1531
1532IF (.NOT.shpisnull(shpobj)) THEN
1533 i = shpwriteobject(shphandle, nshp, shpobj)
1534 CALL shpdestroyobject(shpobj)
1535ENDIF
1536
1537END SUBROUTINE georef_coord_array_export
1538
1539
1550SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1551TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1552CHARACTER(len=*),INTENT(in) :: shpfile
1553
1554REAL(kind=fp_d) :: minb(4), maxb(4)
1555INTEGER :: i, ns, shptype, dbfnf, dbfnr
1556TYPE(shpfileobject) :: shphandle
1557
1558shphandle = shpopen(trim(shpfile), 'rb')
1559IF (shpfileisnull(shphandle)) THEN
1560 ! log here
1561 CALL raise_error()
1562 RETURN
1563ENDIF
1564
1565! get info about file
1566CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1567IF (ns > 0) THEN ! allocate and read the object
1569 DO i = 1, ns
1570 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1571 ENDDO
1572ENDIF
1573
1574CALL shpclose(shphandle)
1575! pack object to save memory
1577
1578END SUBROUTINE arrayof_georef_coord_array_import
1579
1580
1586SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1587TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1588CHARACTER(len=*),INTENT(in) :: shpfile
1589
1590INTEGER :: i
1591TYPE(shpfileobject) :: shphandle
1592
1593IF (this%arraysize > 0) THEN
1594 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1595ELSE
1596 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1597ENDIF
1598IF (shpfileisnull(shphandle)) THEN
1599 ! log here
1600 CALL raise_error()
1601 RETURN
1602ENDIF
1603
1604DO i = 1, this%arraysize
1605 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1606ENDDO
1607
1608CALL shpclose(shphandle)
1609
1610END SUBROUTINE arrayof_georef_coord_array_export
1611#endif
1612
1624FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1625TYPE(georef_coord), INTENT(IN) :: this
1626TYPE(georef_coord_array), INTENT(IN) :: poly
1627LOGICAL :: inside
1628
1629INTEGER :: i
1630
1631inside = .false.
1633IF (.NOT.ALLOCATED(poly%coord)) RETURN
1634! if outside bounding box stop here
1635IF (poly%bbox_updated) THEN
1636 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1637ENDIF
1638
1639IF (ALLOCATED(poly%parts)) THEN
1640 DO i = 1, SIZE(poly%parts)-1
1642 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1643 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1644 ENDDO
1645 IF (SIZE(poly%parts) > 0) THEN ! safety check
1647 poly%coord(poly%parts(i)+1:)%x, &
1648 poly%coord(poly%parts(i)+1:)%y)
1649 ENDIF
1650
1651ELSE
1652 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1653 inside = pointinpoly(this%x, this%y, &
1654 poly%coord(:)%x, poly%coord(:)%y)
1655ENDIF
1656
1657CONTAINS
1658
1659FUNCTION pointinpoly(x, y, px, py)
1660DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1661LOGICAL :: pointinpoly
1662
1663INTEGER :: i, j, starti
1664
1665pointinpoly = .false.
1666
1667IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1668 starti = 2
1669 j = 1
1670ELSE ! unclosed polygon
1671 starti = 1
1672 j = SIZE(px)
1673ENDIF
1674DO i = starti, SIZE(px)
1675 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1676 (py(j) <= y .AND. y < py(i))) THEN
1677 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1678 pointinpoly = .NOT. pointinpoly
1679 ENDIF
1680 ENDIF
1681 j = i
1682ENDDO
1683
1684END FUNCTION pointinpoly
1685
1686END FUNCTION georef_coord_inside
1687
1688
1689
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 |