libsim Versione 7.1.11
|
◆ 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 906 del file georef_coord_class.F90. 907! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
908! authors:
909! Davide Cesari <dcesari@arpa.emr.it>
910! Paolo Patruno <ppatruno@arpa.emr.it>
911
912! This program is free software; you can redistribute it and/or
913! modify it under the terms of the GNU General Public License as
914! published by the Free Software Foundation; either version 2 of
915! the License, or (at your option) any later version.
916
917! This program is distributed in the hope that it will be useful,
918! but WITHOUT ANY WARRANTY; without even the implied warranty of
919! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
920! GNU General Public License for more details.
921
922! You should have received a copy of the GNU General Public License
923! along with this program. If not, see <http://www.gnu.org/licenses/>.
924#include "config.h"
925
943USE geo_proj_class
944#ifdef HAVE_SHAPELIB
945USE shapelib
946#endif
947IMPLICIT NONE
948
954 PRIVATE
955 DOUBLE PRECISION :: x=dmiss, y=dmiss
957
960
966 PRIVATE
967 INTEGER,ALLOCATABLE :: parts(:)
968 TYPE(georef_coord),ALLOCATABLE :: coord(:)
969 INTEGER :: topo=imiss
970 TYPE(geo_proj) :: proj
971 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
972 LOGICAL :: bbox_updated=.false.
974
975INTEGER,PARAMETER :: georef_coord_array_point = 1
976INTEGER,PARAMETER :: georef_coord_array_arc = 3
977INTEGER,PARAMETER :: georef_coord_array_polygon = 5
978INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
979
980
985 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
986END INTERFACE
987
990 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
991END INTERFACE
992
995 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
996END INTERFACE
997
998INTERFACE compute_bbox
999 MODULE PROCEDURE georef_coord_array_compute_bbox
1000END INTERFACE
1001
1003INTERFACE OPERATOR (==)
1004 MODULE PROCEDURE georef_coord_eq
1005END INTERFACE
1006
1008INTERFACE OPERATOR (/=)
1009 MODULE PROCEDURE georef_coord_ne
1010END INTERFACE
1011
1014INTERFACE OPERATOR (>=)
1015 MODULE PROCEDURE georef_coord_ge
1016END INTERFACE
1017
1020INTERFACE OPERATOR (<=)
1021 MODULE PROCEDURE georef_coord_le
1022END INTERFACE
1023
1024#ifdef HAVE_SHAPELIB
1025
1027INTERFACE import
1028 MODULE PROCEDURE arrayof_georef_coord_array_import
1029END INTERFACE
1030
1034 MODULE PROCEDURE arrayof_georef_coord_array_export
1035END INTERFACE
1036#endif
1037
1041 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
1042END INTERFACE
1043
1047 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
1048END INTERFACE
1049
1052 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
1053END INTERFACE
1054
1057 MODULE PROCEDURE georef_coord_dist
1058END INTERFACE
1059
1060#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
1061#define ARRAYOF_TYPE arrayof_georef_coord_array
1062!define ARRAYOF_ORIGEQ 0
1063#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
1064#include "arrayof_pre.F90"
1065! from arrayof
1067
1068PRIVATE
1070 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
1071 georef_coord_array_polygon, georef_coord_array_multipoint, &
1073#ifdef HAVE_SHAPELIB
1075#endif
1077 georef_coord_new, georef_coord_array_new
1078
1079CONTAINS
1080
1081#include "arrayof_post.F90"
1082
1083! ===================
1084! == georef_coord ==
1085! ===================
1089FUNCTION georef_coord_new(x, y) RESULT(this)
1090DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
1091DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
1092TYPE(georef_coord) :: this
1093
1096
1097END FUNCTION georef_coord_new
1098
1099
1100SUBROUTINE georef_coord_delete(this)
1101TYPE(georef_coord),INTENT(inout) :: this
1102
1103this%x = dmiss
1104this%y = dmiss
1105
1106END SUBROUTINE georef_coord_delete
1107
1108
1109ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
1110TYPE(georef_coord),INTENT(in) :: this
1111LOGICAL :: res
1112
1113res = .NOT. this == georef_coord_miss
1114
1115END FUNCTION georef_coord_c_e
1116
1117
1124ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
1125TYPE(georef_coord),INTENT(in) :: this
1126DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1127DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1128
1129IF (PRESENT(x)) x = this%x
1130IF (PRESENT(y)) y = this%y
1131
1132END SUBROUTINE georef_coord_getval
1133
1134
1143ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
1144TYPE(georef_coord),INTENT(in) :: this
1145TYPE(geo_proj),INTENT(in) :: proj
1146DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
1147DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
1148DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
1149DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
1150
1151DOUBLE PRECISION :: llon, llat
1152
1153IF (PRESENT(x)) x = this%x
1154IF (PRESENT(y)) y = this%y
1155IF (PRESENT(lon) .OR. present(lat)) THEN
1157 IF (PRESENT(lon)) lon = llon
1158 IF (PRESENT(lat)) lat = llat
1159ENDIF
1160
1161END SUBROUTINE georef_coord_proj_getval
1162
1163
1164! document and improve
1165ELEMENTAL FUNCTION getlat(this)
1166TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1167DOUBLE PRECISION :: getlat ! latitudine geografica
1168
1169getlat = this%y ! change!!!
1170
1171END FUNCTION getlat
1172
1173! document and improve
1174ELEMENTAL FUNCTION getlon(this)
1175TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
1176DOUBLE PRECISION :: getlon ! longitudine geografica
1177
1178getlon = this%x ! change!!!
1179
1180END FUNCTION getlon
1181
1182
1183ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
1184TYPE(georef_coord),INTENT(IN) :: this, that
1185LOGICAL :: res
1186
1187res = (this%x == that%x .AND. this%y == that%y)
1188
1189END FUNCTION georef_coord_eq
1190
1191
1192ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
1193TYPE(georef_coord),INTENT(IN) :: this, that
1194LOGICAL :: res
1195
1196res = (this%x >= that%x .AND. this%y >= that%y)
1197
1198END FUNCTION georef_coord_ge
1199
1200
1201ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
1202TYPE(georef_coord),INTENT(IN) :: this, that
1203LOGICAL :: res
1204
1205res = (this%x <= that%x .AND. this%y <= that%y)
1206
1207END FUNCTION georef_coord_le
1208
1209
1210ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
1211TYPE(georef_coord),INTENT(IN) :: this, that
1212LOGICAL :: res
1213
1214res = .NOT.(this == that)
1215
1216END FUNCTION georef_coord_ne
1217
1218
1224SUBROUTINE georef_coord_read_unit(this, unit)
1225TYPE(georef_coord),INTENT(out) :: this
1226INTEGER, INTENT(in) :: unit
1227
1228CALL georef_coord_vect_read_unit((/this/), unit)
1229
1230END SUBROUTINE georef_coord_read_unit
1231
1232
1238SUBROUTINE georef_coord_vect_read_unit(this, unit)
1239TYPE(georef_coord) :: this(:)
1240INTEGER, INTENT(in) :: unit
1241
1242CHARACTER(len=40) :: form
1243INTEGER :: i
1244
1245INQUIRE(unit, form=form)
1246IF (form == 'FORMATTED') THEN
1247 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1248!TODO bug gfortran compiler !
1249!missing values are unredeable when formatted
1250ELSE
1251 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1252ENDIF
1253
1254END SUBROUTINE georef_coord_vect_read_unit
1255
1256
1261SUBROUTINE georef_coord_write_unit(this, unit)
1262TYPE(georef_coord),INTENT(in) :: this
1263INTEGER, INTENT(in) :: unit
1264
1265CALL georef_coord_vect_write_unit((/this/), unit)
1266
1267END SUBROUTINE georef_coord_write_unit
1268
1269
1274SUBROUTINE georef_coord_vect_write_unit(this, unit)
1275TYPE(georef_coord),INTENT(in) :: this(:)
1276INTEGER, INTENT(in) :: unit
1277
1278CHARACTER(len=40) :: form
1279INTEGER :: i
1280
1281INQUIRE(unit, form=form)
1282IF (form == 'FORMATTED') THEN
1283 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
1284!TODO bug gfortran compiler !
1285!missing values are unredeable when formatted
1286ELSE
1287 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
1288ENDIF
1289
1290END SUBROUTINE georef_coord_vect_write_unit
1291
1292
1295FUNCTION georef_coord_dist(this, that) RESULT(dist)
1297TYPE(georef_coord), INTENT (IN) :: this
1298TYPE(georef_coord), INTENT (IN) :: that
1299DOUBLE PRECISION :: dist
1300
1301DOUBLE PRECISION :: x,y
1302! Distanza approssimata, valida per piccoli angoli
1303
1304x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
1305y = (this%y-that%y)
1306dist = sqrt(x**2 + y**2)*degrad*rearth
1307
1308END FUNCTION georef_coord_dist
1309
1310
1316FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
1317TYPE(georef_coord),INTENT(IN) :: this
1318TYPE(georef_coord),INTENT(IN) :: coordmin
1319TYPE(georef_coord),INTENT(IN) :: coordmax
1320LOGICAL :: res
1321
1322res = (this >= coordmin .AND. this <= coordmax)
1323
1324END FUNCTION georef_coord_inside_rectang
1325
1326
1327! ========================
1328! == georef_coord_array ==
1329! ========================
1335FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
1336DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
1337DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
1338INTEGER,INTENT(in),OPTIONAL :: topo
1339TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
1340TYPE(georef_coord_array) :: this
1341
1342INTEGER :: lsize
1343
1344IF (PRESENT(x) .AND. PRESENT(y)) THEN
1345 lsize = min(SIZE(x), SIZE(y))
1346 ALLOCATE(this%coord(lsize))
1347 this%coord(1:lsize)%x = x(1:lsize)
1348 this%coord(1:lsize)%y = y(1:lsize)
1349ENDIF
1350this%topo = optio_l(topo)
1352
1353END FUNCTION georef_coord_array_new
1354
1355
1356SUBROUTINE georef_coord_array_delete(this)
1357TYPE(georef_coord_array),INTENT(inout) :: this
1358
1359TYPE(georef_coord_array) :: lobj
1360
1361this = lobj
1362
1363END SUBROUTINE georef_coord_array_delete
1364
1365
1366ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
1367TYPE(georef_coord_array),INTENT(in) :: this
1368LOGICAL :: res
1369
1370res = ALLOCATED(this%coord)
1371
1372END FUNCTION georef_coord_array_c_e
1373
1374
1379SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
1380TYPE(georef_coord_array),INTENT(in) :: this
1381DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
1382DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
1383! allocatable per vedere di nascosto l'effetto che fa
1384INTEGER,OPTIONAL,INTENT(out) :: topo
1385TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
1386
1387
1388IF (PRESENT(x)) THEN
1389 IF (ALLOCATED(this%coord)) THEN
1390 x = this%coord%x
1391 ENDIF
1392ENDIF
1393IF (PRESENT(y)) THEN
1394 IF (ALLOCATED(this%coord)) THEN
1395 y = this%coord%y
1396 ENDIF
1397ENDIF
1398IF (PRESENT(topo)) topo = this%topo
1400
1401END SUBROUTINE georef_coord_array_getval
1402
1403
1409SUBROUTINE georef_coord_array_compute_bbox(this)
1410TYPE(georef_coord_array),INTENT(inout) :: this
1411
1412IF (ALLOCATED(this%coord)) THEN
1413 this%bbox(1)%x = minval(this%coord(:)%x)
1414 this%bbox(1)%y = minval(this%coord(:)%y)
1415 this%bbox(2)%x = maxval(this%coord(:)%x)
1416 this%bbox(2)%y = maxval(this%coord(:)%y)
1417 this%bbox_updated = .true.
1418ENDIF
1419
1420END SUBROUTINE georef_coord_array_compute_bbox
1421
1422#ifdef HAVE_SHAPELIB
1423! internal method for importing a single shape
1424SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
1425TYPE(georef_coord_array),INTENT(OUT) :: this
1426TYPE(shpfileobject),INTENT(INOUT) :: shphandle
1427INTEGER,INTENT(IN) :: nshp
1428
1429TYPE(shpobject) :: shpobj
1430
1431! read shape object
1432shpobj = shpreadobject(shphandle, nshp)
1433IF (.NOT.shpisnull(shpobj)) THEN
1434! import it in georef_coord object
1435 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
1436 topo=shpobj%nshptype)
1437 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
1438 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
1439 ELSE IF (ALLOCATED(this%parts)) THEN
1440 DEALLOCATE(this%parts)
1441 ENDIF
1442 CALL shpdestroyobject(shpobj)
1443 CALL compute_bbox(this)
1444ENDIF
1445
1446
1447END SUBROUTINE georef_coord_array_import
1448
1449
1450! internal method for exporting a single shape
1451SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
1452TYPE(georef_coord_array),INTENT(in) :: this
1453TYPE(shpfileobject),INTENT(inout) :: shphandle
1454INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
1455
1456INTEGER :: i
1457TYPE(shpobject) :: shpobj
1458
1459IF (ALLOCATED(this%coord)) THEN
1460 IF (ALLOCATED(this%parts)) THEN
1461 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
1462 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
1463 ELSE
1464 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
1465 this%coord(:)%x, this%coord(:)%y)
1466 ENDIF
1467ELSE
1468 RETURN
1469ENDIF
1470
1471IF (.NOT.shpisnull(shpobj)) THEN
1472 i = shpwriteobject(shphandle, nshp, shpobj)
1473 CALL shpdestroyobject(shpobj)
1474ENDIF
1475
1476END SUBROUTINE georef_coord_array_export
1477
1478
1489SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
1490TYPE(arrayof_georef_coord_array),INTENT(out) :: this
1491CHARACTER(len=*),INTENT(in) :: shpfile
1492
1493REAL(kind=fp_d) :: minb(4), maxb(4)
1494INTEGER :: i, ns, shptype, dbfnf, dbfnr
1495TYPE(shpfileobject) :: shphandle
1496
1497shphandle = shpopen(trim(shpfile), 'rb')
1498IF (shpfileisnull(shphandle)) THEN
1499 ! log here
1500 CALL raise_error()
1501 RETURN
1502ENDIF
1503
1504! get info about file
1505CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
1506IF (ns > 0) THEN ! allocate and read the object
1508 DO i = 1, ns
1509 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
1510 ENDDO
1511ENDIF
1512
1513CALL shpclose(shphandle)
1514! pack object to save memory
1516
1517END SUBROUTINE arrayof_georef_coord_array_import
1518
1519
1525SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
1526TYPE(arrayof_georef_coord_array),INTENT(in) :: this
1527CHARACTER(len=*),INTENT(in) :: shpfile
1528
1529INTEGER :: i
1530TYPE(shpfileobject) :: shphandle
1531
1532IF (this%arraysize > 0) THEN
1533 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
1534ELSE
1535 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
1536ENDIF
1537IF (shpfileisnull(shphandle)) THEN
1538 ! log here
1539 CALL raise_error()
1540 RETURN
1541ENDIF
1542
1543DO i = 1, this%arraysize
1544 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
1545ENDDO
1546
1547CALL shpclose(shphandle)
1548
1549END SUBROUTINE arrayof_georef_coord_array_export
1550#endif
1551
1563FUNCTION georef_coord_inside(this, poly) RESULT(inside)
1564TYPE(georef_coord), INTENT(IN) :: this
1565TYPE(georef_coord_array), INTENT(IN) :: poly
1566LOGICAL :: inside
1567
1568INTEGER :: i
1569
1570inside = .false.
1572IF (.NOT.ALLOCATED(poly%coord)) RETURN
1573! if outside bounding box stop here
1574IF (poly%bbox_updated) THEN
1575 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
1576ENDIF
1577
1578IF (ALLOCATED(poly%parts)) THEN
1579 DO i = 1, SIZE(poly%parts)-1
1581 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
1582 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
1583 ENDDO
1584 IF (SIZE(poly%parts) > 0) THEN ! safety check
1586 poly%coord(poly%parts(i)+1:)%x, &
1587 poly%coord(poly%parts(i)+1:)%y)
1588 ENDIF
1589
1590ELSE
1591 IF (SIZE(poly%coord) < 1) RETURN ! safety check
1592 inside = pointinpoly(this%x, this%y, &
1593 poly%coord(:)%x, poly%coord(:)%y)
1594ENDIF
1595
1596CONTAINS
1597
1598FUNCTION pointinpoly(x, y, px, py)
1599DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
1600LOGICAL :: pointinpoly
1601
1602INTEGER :: i, j, starti
1603
1604pointinpoly = .false.
1605
1606IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
1607 starti = 2
1608 j = 1
1609ELSE ! unclosed polygon
1610 starti = 1
1611 j = SIZE(px)
1612ENDIF
1613DO i = starti, SIZE(px)
1614 IF ((py(i) <= y .AND. y < py(j)) .OR. &
1615 (py(j) <= y .AND. y < py(i))) THEN
1616 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
1617 pointinpoly = .NOT. pointinpoly
1618 ENDIF
1619 ENDIF
1620 j = i
1621ENDDO
1622
1623END FUNCTION pointinpoly
1624
1625END FUNCTION georef_coord_inside
1626
1627
1628
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 |