libsim Versione 7.2.6
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!> This module defines objects describing georeferenced sparse points
20!! possibly with topology and projection information. This module
21!! defines two classes, \a georef_coord, which represents a single
22!! georeferenced point on the Earth, and \a georef_coord_array which
23!! defines a set of points with a topological relation.
24!!
25!! Both classes have \a PRIVATE members, so that they cannot be
26!! manipulated directly, but only through the proper methods.
27!!
28!! It is also possible to dafine a dynamically extendible array of \a
29!! georef_coord_array objects, of type \a arrayof_georef_coord_array,
30!! suitable for importing/exporting data from/to a shapefile.
31!!
32!! \ingroup base
37USE geo_proj_class
38#ifdef HAVE_SHAPELIB
39USE shapelib
40#endif
41IMPLICIT NONE
42
43!> Derive type defining a single georeferenced point, either in
44!! geodetic or in projected coordinates. The object has no information
45!! about the georeferenced coordinate system associated, it must be
46!! kept separately by the user.
47TYPE georef_coord
48 PRIVATE
49 DOUBLE PRECISION :: x=dmiss, y=dmiss
50END TYPE georef_coord
51
52!> Missing value for georef_coord
53TYPE(georef_coord),PARAMETER :: georef_coord_miss=georef_coord(dmiss,dmiss)
54
55!> Derived type defining a one-dimensional array of georeferenced points
56!! with an associated topology (isolated point, arc, polygon, group of
57!! points), possibly broken into parts and with an associated
58!! georeferenced coordinate system.
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 !< Topology for georef_coord_array (from shapelib): isolated point
70INTEGER,PARAMETER :: georef_coord_array_arc = 3 !< Topology for georef_coord_array (from shapelib): arc (multiple arcs unsupported)
71INTEGER,PARAMETER :: georef_coord_array_polygon = 5 !< Topology for georef_coord_array (from shapelib): polygon (necessarily closed, multiple polygons unsupported)
72INTEGER,PARAMETER :: georef_coord_array_multipoint = 8 !< Topology for georef_coord_array (from shapelib): group of points
73
74
75!> Detructors for the two classes.
76!! They clean up all the information associated with the corresponding
77!! objects.
78INTERFACE delete
79 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
80END INTERFACE
81
82!> Check missing value.
83INTERFACE c_e
84 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
85END INTERFACE
86
87!> Methods for returning the value of object members.
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
96!> Logical equality operator.
97INTERFACE OPERATOR (==)
98 MODULE PROCEDURE georef_coord_eq
99END INTERFACE
100
101!> Logical inequality operator.
102INTERFACE OPERATOR (/=)
103 MODULE PROCEDURE georef_coord_ne
104END INTERFACE
105
106!> Logical greater-equal operator. Returns true if the first point
107!! lies to the northeast of the second
108INTERFACE OPERATOR (>=)
109 MODULE PROCEDURE georef_coord_ge
110END INTERFACE
111
112!> Logical less-equal operator. Returns true if the first point
113!! lies to the southwest of the second
114INTERFACE OPERATOR (<=)
115 MODULE PROCEDURE georef_coord_le
116END INTERFACE
117
118#ifdef HAVE_SHAPELIB
119!> Import an array of \a georef_coord_array objects from a file
120!! in ESRI/Shapefile format.
121INTERFACE import
122 MODULE PROCEDURE arrayof_georef_coord_array_import
123END INTERFACE
124
125!> Export an array of \a georef_coord_array objects to a file
126!! in ESRI/Shapefile format.
127INTERFACE export
128 MODULE PROCEDURE arrayof_georef_coord_array_export
129END INTERFACE
130#endif
131
132!> Read a single \a georef_coord object or an array of \a georef_coord objects
133!! from a Fortran \c FORMATTED or \c UNFORMATTED file.
134INTERFACE read_unit
135 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
136END INTERFACE
137
138!> Write a single \a georef_coord object or an array of \a georef_coord objects
139!! to a Fortran \c FORMATTED or \c UNFORMATTED file.
140INTERFACE write_unit
141 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
142END INTERFACE
143
144!> Determine whether a point lies inside a polygon or a rectangle.
145INTERFACE inside
146 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
147END INTERFACE
148
149!> Compute the distance in m between two points
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! ===================
180!> Construct a \a georef_coord object with the optional parameters provided.
181!! If coordinates are not provided the object obtained is empty
182!! (missing, see c_e function).
183FUNCTION georef_coord_new(x, y) RESULT(this)
184DOUBLE PRECISION,INTENT(in),OPTIONAL :: x !< x coordinate
185DOUBLE PRECISION,INTENT(in),OPTIONAL :: y !< y coordinate
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
212!> Query a \a georef_coord object.
213!! This is the correct way to retrieve the contents of a \a
214!! georef_coord object, since its members are declared as \a
215!! PRIVATE. It is declared as \a ELEMENTAL, thus it works also on
216!! arrays of any shape and, in that case, the result will hae the same
217!! shape as \a this.
218ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
219TYPE(georef_coord),INTENT(in) :: this !< object to query
220DOUBLE PRECISION,INTENT(out),OPTIONAL :: x !< x-coordinate
221DOUBLE PRECISION,INTENT(out),OPTIONAL :: y !< y-coordinate
222
223IF (PRESENT(x)) x = this%x
224IF (PRESENT(y)) y = this%y
225
226END SUBROUTINE georef_coord_getval
227
228
229!> Query a \a georef_coord object associating a geographical projection to it.
230!! This method allow to interpret the \a x,y coordinates of a \a
231!! georef_coord object as projected on a specified geographical
232!! projection system and retrieve the geodetic longitude and/or
233!! latitude associated to them. When \a x or \y are requested it works
234!! as the basic get_val method. It is declared as \a ELEMENTAL, thus
235!! it works also on arrays of any shape and, in that case, the result
236!! will hae the same shape as \a this.
237ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
238TYPE(georef_coord),INTENT(in) :: this !< object to query
239TYPE(geo_proj),INTENT(in) :: proj !< geographical projection associated to coordinates of \a this
240DOUBLE PRECISION,INTENT(out),OPTIONAL :: x !< x-coordinate
241DOUBLE PRECISION,INTENT(out),OPTIONAL :: y !< y-coordinate
242DOUBLE PRECISION,INTENT(out),OPTIONAL :: lon !< geodetic longitude
243DOUBLE PRECISION,INTENT(out),OPTIONAL :: lat !< geodetic latitude
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
313!> Legge da un'unità di file il contenuto dell'oggetto \a this.
314!! Il record da leggere deve essere stato scritto con la ::write_unit
315!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
316!! del vettore devono essere accordate. Il metodo controlla se il file è
317!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
318SUBROUTINE georef_coord_read_unit(this, unit)
319TYPE(georef_coord),INTENT(out) :: this !< oggetto da leggere
320INTEGER, INTENT(in) :: unit !< unità da cui leggere
321
322CALL georef_coord_vect_read_unit((/this/), unit)
323
324END SUBROUTINE georef_coord_read_unit
325
326
327!> Legge da un'unità di file il contenuto dell'oggetto \a this.
328!! Il record da leggere deve essere stato scritto con la ::write_unit
329!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
330!! del vettore devono essere accordate. Il metodo controlla se il file è
331!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
332SUBROUTINE georef_coord_vect_read_unit(this, unit)
333TYPE(georef_coord) :: this(:) !< oggetto da leggere
334INTEGER, INTENT(in) :: unit !< unità da cui leggere
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
351!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
352!! Il record scritto potrà successivamente essere letto con la ::read_unit.
353!! Il metodo controlla se il file è
354!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
355SUBROUTINE georef_coord_write_unit(this, unit)
356TYPE(georef_coord),INTENT(in) :: this !< oggetto da scrivere
357INTEGER, INTENT(in) :: unit !< unità su cui scrivere
358
359CALL georef_coord_vect_write_unit((/this/), unit)
360
361END SUBROUTINE georef_coord_write_unit
362
363
364!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
365!! Il record scritto potrà successivamente essere letto con la ::read_unit.
366!! Il metodo controlla se il file è
367!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
368SUBROUTINE georef_coord_vect_write_unit(this, unit)
369TYPE(georef_coord),INTENT(in) :: this(:) !< oggetto da scrivere
370INTEGER, INTENT(in) :: unit !< unità su cui scrivere
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
387!> Restituisce la distanza in m tra 2 oggetti georef_coord.
388!! La distanza è calcolata approssimativamente ed è valida per piccoli angoli.
389FUNCTION georef_coord_dist(this, that) RESULT(dist)
391TYPE(georef_coord), INTENT (IN) :: this !< primo punto
392TYPE(georef_coord), INTENT (IN) :: that !< secondo punto
393DOUBLE PRECISION :: dist !< distanza in metri
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
405!> Determines whether the point \a this lies inside a specified rectangle.
406!! The rectangle is oriented parallely to the coordinate system, its
407!! lower-left and upper-right vertices are specified by the two
408!! arguments. The function returns \c .TRUE. also in the case it lies
409!! exactly on the border of the rectangle.
410FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
411TYPE(georef_coord),INTENT(IN) :: this !< point to test
412TYPE(georef_coord),INTENT(IN) :: coordmin !< lower-left vertex
413TYPE(georef_coord),INTENT(IN) :: coordmax !< upper-right vertex
414LOGICAL :: res
415
416res = (this >= coordmin .AND. this <= coordmax)
417
418END FUNCTION georef_coord_inside_rectang
419
420
421! ========================
422! == georef_coord_array ==
423! ========================
424!> Construct a \a georef_coord_array object with the optional parameters provided.
425!! If coordinates are not provided the object obtained is empty
426!! (missing, see c_e function), if coordinate arrays are of different
427!! lengths the \a georef_coord_array is initialised to the shortest
428!! length provided.
429FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
430DOUBLE PRECISION,INTENT(in),OPTIONAL :: x(:) !< x coordinate array
431DOUBLE PRECISION,INTENT(in),OPTIONAL :: y(:) !< y coordinate array
432INTEGER,INTENT(in),OPTIONAL :: topo !< topology of the object, one of the \a georef_coord_array_* constants
433TYPE(geo_proj),INTENT(in),OPTIONAL :: proj !< geographical projection associated to the coordinates
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
469!> Query a \a georef_coord_array object.
470!! This is the correct way to retrieve the contents of a \a
471!! georef_coord_array object, since its members are declared as \a
472!! PRIVATE.
473SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
474TYPE(georef_coord_array),INTENT(in) :: this !< object to query
475DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: x(:) !< x-coordinate
476DOUBLE PRECISION,OPTIONAL,ALLOCATABLE,INTENT(out) :: y(:) !< y-coordinate
477! allocatable per vedere di nascosto l'effetto che fa
478INTEGER,OPTIONAL,INTENT(out) :: topo !< topology associated with the coordinates
479TYPE(geo_proj),OPTIONAL,INTENT(out) :: proj !< geographical projection
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
498!> Compute the bounding box of each shape in \a georef_coord_array object.
499!! The bounding box is computed and stored in the object, it is used
500!! by the inside() function for speedup; after it is computed the
501!! object cannot be changed, otherwise the bounding box will not be
502!! valid.
503SUBROUTINE georef_coord_array_compute_bbox(this)
504TYPE(georef_coord_array),INTENT(inout) :: this !< object to manipulate
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
573!> Import an array of \a georef_coord_array objects from a file
574!! in ESRI/Shapefile format. The \a this argument is an uninitialised
575!! \a arrayof_georef_coord_array, every element of which,
576!! this%array(n), is of type \a georef_coord_array and, on return,
577!! will contain information from the n-th shape of the file. Topology
578!! information and possible polygon parts are imported as well, while
579!! no projection information, even if available, is imported. An
580!! error condition while opening the file can be detected by checking
581!! .NOT.ASSOCIATED(this%array), while an error reading shape \a n can
582!! be detected by checking .NOT.c_e(this%array(n)).
583SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
584TYPE(arrayof_georef_coord_array),INTENT(out) :: this !< uninitialised array object
585CHARACTER(len=*),INTENT(in) :: shpfile !< name of shapefile (with or without extension)
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
614!> Export an array of \a georef_coord_array objects to a file
615!! in ESRI/Shapefile format. All the \a this%arraysize shapes
616!! contained in \a this are exported to the requested
617!! shapefile. Topology information and possible polygon parts are
618!! exported as well, while projection information is ignored.
619SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
620TYPE(arrayof_georef_coord_array),INTENT(in) :: this !< array object to be exported
621CHARACTER(len=*),INTENT(in) :: shpfile !< name of shapefile (with or without extension)
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
646!> Determines whether the point \a this lies inside the polygon \a poly.
647!! The polygon is forced to be closed if it is not already the case,
648!! and there is no check about the topology of \a poly to really be of
649!! polygon type. It works also with polygons in parts (as from
650!! shapefile specification) defining either multiple polygons or
651!! polygons with holes.
652!!
653!! The method used consists in counting the number of intersections as
654!! indicated in comp.graphics.algorithms FAQ
655!! (http://www.faqs.org/faqs/graphics/algorithms-faq/) or in
656!! http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html.
657FUNCTION georef_coord_inside(this, poly) RESULT(inside)
658TYPE(georef_coord), INTENT(IN) :: this !< oggetto di cui determinare la posizione
659TYPE(georef_coord_array), INTENT(IN) :: poly !< poligono contenitore
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.