libsim Versione 7.2.1
georef_coord_class.F90
1! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2! authors:
3! Davide Cesari <dcesari@arpa.emr.it>
4! Paolo Patruno <ppatruno@arpa.emr.it>
5
6! This program is free software; you can redistribute it and/or
7! modify it under the terms of the GNU General Public License as
8! published by the Free Software Foundation; either version 2 of
9! the License, or (at your option) any later version.
10
11! This program is distributed in the hope that it will be useful,
12! but WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14! GNU General Public License for more details.
15
16! You should have received a copy of the GNU General Public License
17! along with this program. If not, see <http://www.gnu.org/licenses/>.
18#include "config.h"
19
37USE geo_proj_class
38#ifdef HAVE_SHAPELIB
39USE shapelib
40#endif
41IMPLICIT NONE
42
47TYPE georef_coord
48 PRIVATE
49 DOUBLE PRECISION :: x=dmiss, y=dmiss
50END TYPE georef_coord
51
53TYPE(georef_coord),PARAMETER :: georef_coord_miss=georef_coord(dmiss,dmiss)
54
60 PRIVATE
61 INTEGER,ALLOCATABLE :: parts(:)
62 TYPE(georef_coord),ALLOCATABLE :: coord(:)
63 INTEGER :: topo=imiss
64 TYPE(geo_proj) :: proj
65 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
66 LOGICAL :: bbox_updated=.false.
67END TYPE georef_coord_array
68
69INTEGER,PARAMETER :: georef_coord_array_point = 1
70INTEGER,PARAMETER :: georef_coord_array_arc = 3
71INTEGER,PARAMETER :: georef_coord_array_polygon = 5
72INTEGER,PARAMETER :: georef_coord_array_multipoint = 8
73
74
78INTERFACE delete
79 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
80END INTERFACE
81
83INTERFACE c_e
84 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
85END INTERFACE
86
88INTERFACE getval
89 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
90END INTERFACE
91
92INTERFACE compute_bbox
93 MODULE PROCEDURE georef_coord_array_compute_bbox
94END INTERFACE
95
97INTERFACE OPERATOR (==)
98 MODULE PROCEDURE georef_coord_eq
99END INTERFACE
100
102INTERFACE OPERATOR (/=)
103 MODULE PROCEDURE georef_coord_ne
104END INTERFACE
105
108INTERFACE OPERATOR (>=)
109 MODULE PROCEDURE georef_coord_ge
110END INTERFACE
111
114INTERFACE OPERATOR (<=)
115 MODULE PROCEDURE georef_coord_le
116END INTERFACE
117
118#ifdef HAVE_SHAPELIB
119
121INTERFACE import
122 MODULE PROCEDURE arrayof_georef_coord_array_import
123END INTERFACE
124
127INTERFACE export
128 MODULE PROCEDURE arrayof_georef_coord_array_export
129END INTERFACE
130#endif
131
134INTERFACE read_unit
135 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
136END INTERFACE
137
140INTERFACE write_unit
141 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
142END INTERFACE
143
145INTERFACE inside
146 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
147END INTERFACE
148
150INTERFACE dist
151 MODULE PROCEDURE georef_coord_dist
152END INTERFACE
153
154#define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
155#define ARRAYOF_TYPE arrayof_georef_coord_array
156!define ARRAYOF_ORIGEQ 0
157#define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
158#include "arrayof_pre.F90"
159! from arrayof
161
162PRIVATE
163PUBLIC georef_coord, georef_coord_miss, &
164 georef_coord_array, georef_coord_array_point, georef_coord_array_arc, &
165 georef_coord_array_polygon, georef_coord_array_multipoint, &
166 delete, c_e, getval, compute_bbox, OPERATOR(==), OPERATOR(/=), OPERATOR(>=), OPERATOR(<=), &
167#ifdef HAVE_SHAPELIB
168 import, export, &
169#endif
171 georef_coord_new, georef_coord_array_new
172
173CONTAINS
174
175#include "arrayof_post.F90"
176
177! ===================
178! == georef_coord ==
179! ===================
183FUNCTION georef_coord_new(x, y) RESULT(this)
184DOUBLE PRECISION,INTENT(in),OPTIONAL :: x
185DOUBLE PRECISION,INTENT(in),OPTIONAL :: y
186TYPE(georef_coord) :: this
187
188CALL optio(x, this%x)
189CALL optio(y, this%y)
190
191END FUNCTION georef_coord_new
192
193
194SUBROUTINE georef_coord_delete(this)
195TYPE(georef_coord),INTENT(inout) :: this
196
197this%x = dmiss
198this%y = dmiss
199
200END SUBROUTINE georef_coord_delete
201
202
203ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
204TYPE(georef_coord),INTENT(in) :: this
205LOGICAL :: res
206
207res = .NOT. this == georef_coord_miss
208
209END FUNCTION georef_coord_c_e
210
211
218ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
219TYPE(georef_coord),INTENT(in) :: this
220DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
221DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
222
223IF (PRESENT(x)) x = this%x
224IF (PRESENT(y)) y = this%y
225
226END SUBROUTINE georef_coord_getval
227
228
237ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
238TYPE(georef_coord),INTENT(in) :: this
239TYPE(geo_proj),INTENT(in) :: proj
240DOUBLE PRECISION,INTENT(out),OPTIONAL :: x
241DOUBLE PRECISION,INTENT(out),OPTIONAL :: y
242DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon
243DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat
244
245DOUBLE PRECISION :: llon, llat
246
247IF (PRESENT(x)) x = this%x
248IF (PRESENT(y)) y = this%y
249IF (PRESENT(lon) .OR. present(lat)) THEN
250 CALL unproj(proj, this%x, this%y, llon, llat)
251 IF (PRESENT(lon)) lon = llon
252 IF (PRESENT(lat)) lat = llat
253ENDIF
254
255END SUBROUTINE georef_coord_proj_getval
256
258! document and improve
259ELEMENTAL FUNCTION getlat(this)
260TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
261DOUBLE PRECISION :: getlat ! latitudine geografica
262
263getlat = this%y ! change!!!
264
265END FUNCTION getlat
267! document and improve
268ELEMENTAL FUNCTION getlon(this)
269TYPE(georef_coord),INTENT(in) :: this ! oggetto di cui restituire latitudine
270DOUBLE PRECISION :: getlon ! longitudine geografica
272getlon = this%x ! change!!!
273
274END FUNCTION getlon
275
277ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
278TYPE(georef_coord),INTENT(IN) :: this, that
279LOGICAL :: res
280
281res = (this%x == that%x .AND. this%y == that%y)
282
283END FUNCTION georef_coord_eq
284
286ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
287TYPE(georef_coord),INTENT(IN) :: this, that
288LOGICAL :: res
289
290res = (this%x >= that%x .AND. this%y >= that%y)
291
292END FUNCTION georef_coord_ge
293
294
295ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
296TYPE(georef_coord),INTENT(IN) :: this, that
297LOGICAL :: res
298
299res = (this%x <= that%x .AND. this%y <= that%y)
300
301END FUNCTION georef_coord_le
303
304ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
305TYPE(georef_coord),INTENT(IN) :: this, that
306LOGICAL :: res
307
308res = .NOT.(this == that)
310END FUNCTION georef_coord_ne
311
312
318SUBROUTINE georef_coord_read_unit(this, unit)
319TYPE(georef_coord),INTENT(out) :: this
320INTEGER, INTENT(in) :: unit
321
322CALL georef_coord_vect_read_unit((/this/), unit)
323
324END SUBROUTINE georef_coord_read_unit
325
326
332SUBROUTINE georef_coord_vect_read_unit(this, unit)
333TYPE(georef_coord) :: this(:)
334INTEGER, INTENT(in) :: unit
335
336CHARACTER(len=40) :: form
337INTEGER :: i
339INQUIRE(unit, form=form)
340IF (form == 'FORMATTED') THEN
341 read(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
342!TODO bug gfortran compiler !
343!missing values are unredeable when formatted
344ELSE
345 READ(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
346ENDIF
347
348END SUBROUTINE georef_coord_vect_read_unit
349
350
355SUBROUTINE georef_coord_write_unit(this, unit)
356TYPE(georef_coord),INTENT(in) :: this
357INTEGER, INTENT(in) :: unit
358
359CALL georef_coord_vect_write_unit((/this/), unit)
360
361END SUBROUTINE georef_coord_write_unit
362
363
368SUBROUTINE georef_coord_vect_write_unit(this, unit)
369TYPE(georef_coord),INTENT(in) :: this(:)
370INTEGER, INTENT(in) :: unit
371
372CHARACTER(len=40) :: form
373INTEGER :: i
374
375INQUIRE(unit, form=form)
376IF (form == 'FORMATTED') THEN
377 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,SIZE(this))
378!TODO bug gfortran compiler !
379!missing values are unredeable when formatted
380ELSE
381 WRITE(unit) (this(i)%x,this(i)%y, i=1,SIZE(this))
382ENDIF
383
384END SUBROUTINE georef_coord_vect_write_unit
385
389FUNCTION georef_coord_dist(this, that) RESULT(dist)
391TYPE(georef_coord), INTENT (IN) :: this
392TYPE(georef_coord), INTENT (IN) :: that
393DOUBLE PRECISION :: dist
394
395DOUBLE PRECISION :: x,y
396! Distanza approssimata, valida per piccoli angoli
397
398x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
399y = (this%y-that%y)
400dist = sqrt(x**2 + y**2)*degrad*rearth
402END FUNCTION georef_coord_dist
403
404
410FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
411TYPE(georef_coord),INTENT(IN) :: this
412TYPE(georef_coord),INTENT(IN) :: coordmin
413TYPE(georef_coord),INTENT(IN) :: coordmax
414LOGICAL :: res
415
416res = (this >= coordmin .AND. this <= coordmax)
417
418END FUNCTION georef_coord_inside_rectang
419
420
421! ========================
422! == georef_coord_array ==
423! ========================
429FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
430DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:)
431DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:)
432INTEGER,INTENT(in),OPTIONAL :: topo
433TYPE(geo_proj),INTENT(in),OPTIONAL :: proj
434TYPE(georef_coord_array) :: this
435
436INTEGER :: lsize
437
438IF (PRESENT(x) .AND. PRESENT(y)) THEN
439 lsize = min(SIZE(x), SIZE(y))
440 ALLOCATE(this%coord(lsize))
441 this%coord(1:lsize)%x = x(1:lsize)
442 this%coord(1:lsize)%y = y(1:lsize)
443ENDIF
444this%topo = optio_l(topo)
445IF (PRESENT(proj)) this%proj = proj
446
447END FUNCTION georef_coord_array_new
448
449
450SUBROUTINE georef_coord_array_delete(this)
451TYPE(georef_coord_array),INTENT(inout) :: this
452
453TYPE(georef_coord_array) :: lobj
454
455this = lobj
456
457END SUBROUTINE georef_coord_array_delete
458
459
460ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
461TYPE(georef_coord_array),INTENT(in) :: this
462LOGICAL :: res
463
464res = ALLOCATED(this%coord)
465
466END FUNCTION georef_coord_array_c_e
467
468
473SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
474TYPE(georef_coord_array),INTENT(in) :: this
475DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:)
476DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:)
477! allocatable per vedere di nascosto l'effetto che fa
478INTEGER,OPTIONAL,INTENT(out) :: topo
479TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj
480
481
482IF (PRESENT(x)) THEN
483 IF (ALLOCATED(this%coord)) THEN
484 x = this%coord%x
485 ENDIF
486ENDIF
487IF (PRESENT(y)) THEN
488 IF (ALLOCATED(this%coord)) THEN
489 y = this%coord%y
490 ENDIF
491ENDIF
492IF (PRESENT(topo)) topo = this%topo
493IF (PRESENT(proj)) proj = this%proj ! warning proj has no missing value yet
494
495END SUBROUTINE georef_coord_array_getval
496
497
503SUBROUTINE georef_coord_array_compute_bbox(this)
504TYPE(georef_coord_array),INTENT(inout) :: this
505
506IF (ALLOCATED(this%coord)) THEN
507 this%bbox(1)%x = minval(this%coord(:)%x)
508 this%bbox(1)%y = minval(this%coord(:)%y)
509 this%bbox(2)%x = maxval(this%coord(:)%x)
510 this%bbox(2)%y = maxval(this%coord(:)%y)
511 this%bbox_updated = .true.
512ENDIF
513
514END SUBROUTINE georef_coord_array_compute_bbox
516#ifdef HAVE_SHAPELIB
517! internal method for importing a single shape
518SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
519TYPE(georef_coord_array),INTENT(OUT) :: this
520TYPE(shpfileobject),INTENT(INOUT) :: shphandle
521INTEGER,INTENT(IN) :: nshp
522
523TYPE(shpobject) :: shpobj
524
525! read shape object
526shpobj = shpreadobject(shphandle, nshp)
527IF (.NOT.shpisnull(shpobj)) THEN
528! import it in georef_coord object
529 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
530 topo=shpobj%nshptype)
531 IF (shpobj%nparts > 1 .AND. ASSOCIATED(shpobj%panpartstart)) THEN
532 this%parts = shpobj%panpartstart(:) ! automatic f95 allocation
533 ELSE IF (ALLOCATED(this%parts)) THEN
534 DEALLOCATE(this%parts)
535 ENDIF
536 CALL shpdestroyobject(shpobj)
537 CALL compute_bbox(this)
538ENDIF
539
540
541END SUBROUTINE georef_coord_array_import
542
543
544! internal method for exporting a single shape
545SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
546TYPE(georef_coord_array),INTENT(in) :: this
547TYPE(shpfileobject),INTENT(inout) :: shphandle
548INTEGER,INTENT(IN) :: nshp ! index of shape to write starting from 0, -1 to append
549
550INTEGER :: i
551TYPE(shpobject) :: shpobj
552
553IF (ALLOCATED(this%coord)) THEN
554 IF (ALLOCATED(this%parts)) THEN
555 shpobj = shpcreateobject(this%topo, -1, SIZE(this%parts), this%parts, &
556 this%parts, SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
557 ELSE
558 shpobj = shpcreatesimpleobject(this%topo, SIZE(this%coord), &
559 this%coord(:)%x, this%coord(:)%y)
560 ENDIF
561ELSE
562 RETURN
563ENDIF
564
565IF (.NOT.shpisnull(shpobj)) THEN
566 i = shpwriteobject(shphandle, nshp, shpobj)
567 CALL shpdestroyobject(shpobj)
568ENDIF
569
570END SUBROUTINE georef_coord_array_export
571
572
583SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
584TYPE(arrayof_georef_coord_array),INTENT(out) :: this
585CHARACTER(len=*),INTENT(in) :: shpfile
586
587REAL(kind=fp_d) :: minb(4), maxb(4)
588INTEGER :: i, ns, shptype, dbfnf, dbfnr
589TYPE(shpfileobject) :: shphandle
590
591shphandle = shpopen(trim(shpfile), 'rb')
592IF (shpfileisnull(shphandle)) THEN
593 ! log here
594 CALL raise_error()
595 RETURN
596ENDIF
597
598! get info about file
599CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
600IF (ns > 0) THEN ! allocate and read the object
601 CALL insert(this, nelem=ns)
602 DO i = 1, ns
603 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
604 ENDDO
605ENDIF
606
607CALL shpclose(shphandle)
608! pack object to save memory
609CALL packarray(this)
610
611END SUBROUTINE arrayof_georef_coord_array_import
612
613
619SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
620TYPE(arrayof_georef_coord_array),INTENT(in) :: this
621CHARACTER(len=*),INTENT(in) :: shpfile
622
623INTEGER :: i
624TYPE(shpfileobject) :: shphandle
625
626IF (this%arraysize > 0) THEN
627 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
628ELSE
629 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
630ENDIF
631IF (shpfileisnull(shphandle)) THEN
632 ! log here
633 CALL raise_error()
634 RETURN
635ENDIF
636
637DO i = 1, this%arraysize
638 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
639ENDDO
640
641CALL shpclose(shphandle)
642
643END SUBROUTINE arrayof_georef_coord_array_export
644#endif
645
657FUNCTION georef_coord_inside(this, poly) RESULT(inside)
658TYPE(georef_coord), INTENT(IN) :: this
659TYPE(georef_coord_array), INTENT(IN) :: poly
660LOGICAL :: inside
661
662INTEGER :: i
663
664inside = .false.
665IF (.NOT.c_e(this)) RETURN
666IF (.NOT.ALLOCATED(poly%coord)) RETURN
667! if outside bounding box stop here
668IF (poly%bbox_updated) THEN
669 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2))) RETURN
670ENDIF
671
672IF (ALLOCATED(poly%parts)) THEN
673 DO i = 1, SIZE(poly%parts)-1
674 inside = inside .NEQV. pointinpoly(this%x, this%y, &
675 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
676 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
677 ENDDO
678 IF (SIZE(poly%parts) > 0) THEN ! safety check
679 inside = inside .NEQV. pointinpoly(this%x, this%y, &
680 poly%coord(poly%parts(i)+1:)%x, &
681 poly%coord(poly%parts(i)+1:)%y)
682 ENDIF
683
684ELSE
685 IF (SIZE(poly%coord) < 1) RETURN ! safety check
686 inside = pointinpoly(this%x, this%y, &
687 poly%coord(:)%x, poly%coord(:)%y)
688ENDIF
689
690CONTAINS
691
692FUNCTION pointinpoly(x, y, px, py)
693DOUBLE PRECISION, INTENT(in) :: x, y, px(:), py(:)
694LOGICAL :: pointinpoly
695
696INTEGER :: i, j, starti
697
698pointinpoly = .false.
699
700IF (px(1) == px(SIZE(px)) .AND. py(1) == py(SIZE(px))) THEN ! closed polygon
701 starti = 2
702 j = 1
703ELSE ! unclosed polygon
704 starti = 1
705 j = SIZE(px)
706ENDIF
707DO i = starti, SIZE(px)
708 IF ((py(i) <= y .AND. y < py(j)) .OR. &
709 (py(j) <= y .AND. y < py(i))) THEN
710 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i)) THEN
711 pointinpoly = .NOT. pointinpoly
712 ENDIF
713 ENDIF
714 j = i
715ENDDO
716
717END FUNCTION pointinpoly
718
719END FUNCTION georef_coord_inside
720
721
722
723END MODULE georef_coord_class
Compute forward coordinate transformation from geographical system to projected system.
Compute backward coordinate transformation from projected system to geographical system.
Quick method to append an element to the array.
Detructors for the two classes.
Compute the distance in m between two points.
Export an array of georef_coord_array objects to a file in ESRI/Shapefile format.
Methods for returning the value of object members.
Import an array of georef_coord_array objects from a file in ESRI/Shapefile format.
Method for inserting elements of the array at a desired position.
Determine whether a point lies inside a polygon or a rectangle.
Method for packing the array object reducing at a minimum the memory occupation, without destroying i...
Read a single georef_coord object or an array of georef_coord objects from a Fortran FORMATTED or UNF...
Method for removing elements of the array at a desired position.
Write a single georef_coord object or an array of georef_coord objects to a Fortran FORMATTED or UNFO...
Generic subroutine for checking OPTIONAL parameters.
Costanti fisiche (DOUBLEPRECISION).
Gestione degli errori.
This module defines objects describing georeferenced sparse points possibly with topology and project...
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Derived type defining a dynamically extensible array of TYPE(georef_coord_array) elements.
Derived type defining a one-dimensional array of georeferenced points with an associated topology (is...
Derive type defining a single georeferenced point, either in geodetic or in projected coordinates.

Generated with Doxygen.