libsim Versione 7.2.6
geo_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!> Classes for handling georeferenced sparse points in geographical
20!! corodinates. This module defines two classes for managing
21!! georeferenced points on the Earth in geographical polar
22!! coordinates. It allows importing and exporting blocks of points
23!! from/to plain text files and \a ESRI \a shapefile's. Both classes
24!! have \a PRIVATE members, so that they cannot be manipulated
25!! directly, but only through the proper methods.
26!!
27!! \ingroup base
28MODULE geo_coord_class
29USE kinds
33!USE doubleprecision_phys_const
35#ifdef HAVE_SHAPELIB
36USE shapelib
37!, ONLY: shpobject, shpreadobject, shpdestroyobject, shpwriteobject, &
38! shpcreatesimpleobject
39#endif
40IMPLICIT NONE
41
42
43!> REAL Kind for geographical coordinates.
44!! This constant has to be used when defining or converting values to be
45!! used as geographical coordinates, e.g.:
46!! \code
47!! REAL(kind=fp_geo) :: x, y
48!! coordx = REAL(mylon, kind=fp_geo)
49!! \endcode
50INTEGER, PARAMETER :: fp_geo=fp_d
51real(kind=fp_geo), PARAMETER :: fp_geo_miss=dmiss
52
53!> Derived type defining an isolated georeferenced point on Earth in
54!! polar geographical coordinates
55TYPE geo_coord
56 PRIVATE
57 INTEGER(kind=int_l) :: ilon, ilat
58END TYPE geo_coord
59
60TYPE geo_coord_r
61 PRIVATE
62 REAL(kind=fp_geo) :: lon, lat
63END TYPE geo_coord_r
64
65
66!> Derived type defining a one-dimensional array of georeferenced points
67!! with an associated topology (isolated point, arc, polygon, group of
68!! points)
70 PRIVATE
71 REAL(kind=fp_geo),POINTER :: ll(:,:) => null()
72 INTEGER :: vsize, vtype
73END TYPE geo_coordvect
74
75!> Missing value for geo_coord
76TYPE(geo_coord),PARAMETER :: geo_coord_miss=geo_coord(imiss,imiss)
77
78INTEGER, PARAMETER :: geo_coordvect_point = 1 !< Topology for geo_coordvect (from shapelib): isolated point
79INTEGER, PARAMETER :: geo_coordvect_arc = 3 !< Topology for geo_coordvect (from shapelib): arc (multiple arcs unsupported)
80INTEGER, PARAMETER :: geo_coordvect_polygon = 5 !< Topology for geo_coordvect (from shapelib): polygon (necessarily closed, multiple polygons unsupported)
81INTEGER, PARAMETER :: geo_coordvect_multipoint = 8 !< Topology for geo_coordvect (from shapelib): group of points
82
83
84!> Constructors for the two classes.
85!! They have to be called for every object of these types which is
86!! going to be used.
87INTERFACE init
88 MODULE PROCEDURE geo_coord_init, geo_coordvect_init
89END INTERFACE
90
91!> Detructors for the two classes.
92!! They clean up all the information associated with the corresponding
93!! objects.
94INTERFACE delete
95 MODULE PROCEDURE geo_coord_delete, geo_coordvect_delete
96END INTERFACE
97
98!> Methods for returning the value of object members.
99INTERFACE getval
100 MODULE PROCEDURE geo_coord_getval, geo_coordvect_getval
101END INTERFACE
102
103!> Logical equality operator.
104INTERFACE OPERATOR (==)
105 MODULE PROCEDURE geo_coord_eq
106END INTERFACE
107
108!> Logical inequality operator.
109INTERFACE OPERATOR (/=)
110 MODULE PROCEDURE geo_coord_ne
111END INTERFACE
112
113INTERFACE OPERATOR (>=)
114 MODULE PROCEDURE geo_coord_ge
115END INTERFACE
116
117INTERFACE OPERATOR (>)
118 MODULE PROCEDURE geo_coord_gt
119END INTERFACE
120
121INTERFACE OPERATOR (<=)
122 MODULE PROCEDURE geo_coord_le
123END INTERFACE
124
125INTERFACE OPERATOR (<)
126 MODULE PROCEDURE geo_coord_lt
127END INTERFACE
128
129!> Import one or more \a geo_coordvect objects from a plain text file
130!! or for a file in ESRI/Shapefile format.
131INTERFACE import
132 MODULE PROCEDURE geo_coordvect_import, geo_coordvect_importvect
133END INTERFACE
134
135!> Export one or more \a geo_coordvect objects to a plain text file
136!! or to a file in ESRI/Shapefile format.
137INTERFACE export
138 MODULE PROCEDURE geo_coordvect_export, geo_coordvect_exportvect
139END INTERFACE
140
141!> Read a single \a geo_coord object or an array of \a geo_coord objects
142!! from a Fortran \c FORMATTED or \c UNFORMATTED file.
143INTERFACE read_unit
144 MODULE PROCEDURE geo_coord_read_unit, geo_coord_vect_read_unit
145END INTERFACE
146
147!> Write a single \a geo_coord object or an array of \a geo_coord objects
148!! to a Fortran \c FORMATTED or \c UNFORMATTED file.
149INTERFACE write_unit
150 MODULE PROCEDURE geo_coord_write_unit, geo_coord_vect_write_unit
151END INTERFACE
152
153!> Determine whether a point lies inside a polygon or a rectangle.
154INTERFACE inside
155 MODULE PROCEDURE geo_coord_inside, geo_coord_inside_rectang
156END INTERFACE
157
158!> Missing check
159INTERFACE c_e
160 MODULE PROCEDURE c_e_geo_coord
161END INTERFACE
162
163!>Represent geo_coord object in a pretty string
164INTERFACE to_char
165 MODULE PROCEDURE to_char_geo_coord
166END INTERFACE
167
168!>Print object
169INTERFACE display
170 MODULE PROCEDURE display_geo_coord
171END INTERFACE
172
173CONTAINS
174
175
176! ===================
177! == geo_coord ==
178! ===================
179!> Costruisce un oggetto \a geo_coord con i parametri opzionali forniti.
180!! Se sono presenti \a lon e \a lat, inizializza le coordinate geografiche
181!! ignorando \a utme e \a utmn, mentre se sono specificati \a utme e \a utmn
182!! succede il contrario; non è possibile specificare le coordinate in entrambi
183!! i sistemi, usare eventualmente \a to_geo. Se non viene passato nessun parametro
184!! opzionale l'oggetto è inizializzato a valore mancante.
185SUBROUTINE geo_coord_init(this, lon, lat, ilon, ilat)
186TYPE(geo_coord) :: this !< oggetto da inizializzare
187REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon !< longitudine geografica
188REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat !< latitudine geografica
189integer(kind=int_l), INTENT(IN), OPTIONAL :: ilon !< integer longitudine geografica (nint(lon*1.e5)
190integer(kind=int_l), INTENT(IN), OPTIONAL :: ilat !< integer latitudine geografica (nint(lat*1.e5)
191
192real(kind=fp_geo) :: llon,llat
193
194CALL optio(ilon, this%ilon)
195CALL optio(ilat, this%ilat)
196
197if (.not. c_e(this%ilon)) then
198 CALL optio(lon, llon)
199 if (c_e(llon)) this%ilon=nint(llon*1.d5)
200end if
201
202if (.not. c_e(this%ilat)) then
203 CALL optio(lat, llat)
204 if (c_e(llat)) this%ilat=nint(llat*1.d5)
205end if
206
207END SUBROUTINE geo_coord_init
208
209!> Distrugge l'oggetto in maniera pulita, assegnandogli un valore mancante.
210SUBROUTINE geo_coord_delete(this)
211TYPE(geo_coord), INTENT(INOUT) :: this !< oggetto da distruggre
212
213this%ilon = imiss
214this%ilat = imiss
215
216END SUBROUTINE geo_coord_delete
217
218!> Restituisce il valore di uno o più componenti di un oggetto \a geo_coord.
219!! Qualsiasi combinazione dei parametri opzionali è consentita; se
220!! il tipo di coordinata richiesta non è stato inizializzato né calcolato,
221!! restituisce il corrispondente valore mancante.
222elemental SUBROUTINE geo_coord_getval(this, lon, lat,ilon,ilat)
223TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire i componenti
224REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lon !< longitudine geografica
225REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lat !< latitudine geografica
226integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilon !< integer longitudine geografica (nint(lon*1.e5)
227integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilat !< integer latitudine geografica (nint(lat*1.e5)
228
229IF (PRESENT(ilon)) ilon = getilon(this)
230IF (PRESENT(ilat)) ilat = getilat(this)
231
232IF (PRESENT(lon)) lon = getlon(this)
233IF (PRESENT(lat)) lat = getlat(this)
234
235END SUBROUTINE geo_coord_getval
236
237
238!> Restituisce la latitudine di uno o più componenti di un oggetto \a geo_coord.
239!! Se la latitudine non è stata inizializzata né calcolata
240!! restituisce il corrispondente valore mancante.
241!! Nata per permettere operazioni vettorizzate
242elemental FUNCTION getilat(this)
243TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire latitudine
244integer(kind=int_l) :: getilat !< latitudine geografica
245
246getilat = this%ilat
247
248END FUNCTION getilat
249
250!> Restituisce la latitudine di uno o più componenti di un oggetto \a geo_coord.
251!! Se la latitudine non è stata inizializzata né calcolata
252!! restituisce il corrispondente valore mancante.
253!! Nata per permettere operazioni vettorizzate
254elemental FUNCTION getlat(this)
255TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire latitudine
256real(kind=fp_geo) :: getlat !< latitudine geografica
257integer(kind=int_l) :: ilat
258
259ilat=getilat(this)
260if (c_e(ilat)) then
261 getlat = ilat*1.d-5
262else
263 getlat=fp_geo_miss
264end if
265
266END FUNCTION getlat
268!> Restituisce la longitudine di uno o più componenti di un oggetto \a geo_coord.
269!! Se la latitudine non è stata inizializzata né calcolata
270!! restituisce il corrispondente valore mancante.
271!! Nata per permettere operazioni vettorizzate
272elemental FUNCTION getilon(this)
273TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire latitudine
274integer(kind=int_l) :: getilon !< longitudine geografica
276getilon = this%ilon
277
278END FUNCTION getilon
279
280
281!> Restituisce la longitudine di uno o più componenti di un oggetto \a geo_coord.
282!! Se la latitudine non è stata inizializzata né calcolata
283!! restituisce il corrispondente valore mancante.
284!! Nata per permettere operazioni vettorizzate
285elemental FUNCTION getlon(this)
286TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui restituire latitudine
287real(kind=fp_geo) :: getlon !< longitudine geografica
288integer(kind=int_l) :: ilon
289
290ilon=getilon(this)
291if (c_e(ilon)) then
292 getlon = ilon*1.d-5
293else
294 getlon=fp_geo_miss
295end if
296
297END FUNCTION getlon
298
299
300elemental FUNCTION geo_coord_eq(this, that) RESULT(res)
301TYPE(geo_coord),INTENT(IN) :: this, that
302LOGICAL :: res
303
304res = (this%ilon == that%ilon .AND. this%ilat == that%ilat)
305
306END FUNCTION geo_coord_eq
307
308
309elemental FUNCTION geo_coord_ne(this, that) RESULT(res)
310TYPE(geo_coord),INTENT(IN) :: this, that
311LOGICAL :: res
312
313res = .not. this == that
314
315END FUNCTION geo_coord_ne
316
317!> Logical great operator. Returns true if the first point
318!! lies to the west of the second or if lon is equal north
319elemental FUNCTION geo_coord_gt(this, that) RESULT(res)
320TYPE(geo_coord),INTENT(IN) :: this, that
321LOGICAL :: res
322
323res = this%ilon > that%ilon
324
325if ( this%ilon == that%ilon ) then
326 res = this%ilat > that%ilat
327end if
328
329END FUNCTION geo_coord_gt
330
332!> Logical great-equal operator. Returns true if the first point
333!! lies to the west of the second or if lon is equal north
334elemental FUNCTION geo_coord_ge(this, that) RESULT(res)
335TYPE(geo_coord),INTENT(IN) :: this, that
336LOGICAL :: res
338res = .not. this < that
339
340END FUNCTION geo_coord_ge
341
342!> Logical less operator. Returns true if the first point
343!! lies to the west of the second or if lon is equal south
344elemental FUNCTION geo_coord_lt(this, that) RESULT(res)
345TYPE(geo_coord),INTENT(IN) :: this, that
346LOGICAL :: res
348res = this%ilon < that%ilon
349
350if ( this%ilon == that%ilon ) then
351 res = this%ilat < that%ilat
352end if
353
354
355END FUNCTION geo_coord_lt
356
358!> Logical less-equal operator. Returns true if the first point
359!! lies to the west of the second or if lon is equal south
360elemental FUNCTION geo_coord_le(this, that) RESULT(res)
361TYPE(geo_coord),INTENT(IN) :: this, that
362LOGICAL :: res
363
364res = .not. this > that
365
366END FUNCTION geo_coord_le
367
368
369!> Logical greater-equal operator. Returns true if the first point
370!! lies to the northeast of the second
371elemental FUNCTION geo_coord_ure(this, that) RESULT(res)
372TYPE(geo_coord),INTENT(IN) :: this, that
373LOGICAL :: res
374
375res = (this%ilon >= that%ilon .AND. this%ilat >= that%ilat)
376
377END FUNCTION geo_coord_ure
378
379!> Logical greater operator. Returns true if the first point
380!! lies to the northeast of the second
381elemental FUNCTION geo_coord_ur(this, that) RESULT(res)
382TYPE(geo_coord),INTENT(IN) :: this, that
383LOGICAL :: res
384
385res = (this%ilon > that%ilon .AND. this%ilat > that%ilat)
386
387END FUNCTION geo_coord_ur
388
389
390!> Logical less-equal operator. Returns true if the first point
391!! lies to the southwest of the second
392elemental FUNCTION geo_coord_lle(this, that) RESULT(res)
393TYPE(geo_coord),INTENT(IN) :: this, that
394LOGICAL :: res
395
396res = (this%ilon <= that%ilon .AND. this%ilat <= that%ilat)
397
398END FUNCTION geo_coord_lle
399
400!> Logical less operator. Returns true if the first point
401!! lies to the southwest of the second
402elemental FUNCTION geo_coord_ll(this, that) RESULT(res)
403TYPE(geo_coord),INTENT(IN) :: this, that
404LOGICAL :: res
405
406res = (this%ilon < that%ilon .AND. this%ilat < that%ilat)
407
408END FUNCTION geo_coord_ll
409
411!> Legge da un'unità di file il contenuto dell'oggetto \a this.
412!! Il record da leggere deve essere stato scritto con la ::write_unit
413!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
414!! del vettore devono essere accordate. Il metodo controlla se il file è
415!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
416SUBROUTINE geo_coord_read_unit(this, unit)
417TYPE(geo_coord),INTENT(out) :: this !< oggetto da leggere
418INTEGER, INTENT(in) :: unit !< unità da cui leggere
419
420CALL geo_coord_vect_read_unit((/this/), unit)
421
422END SUBROUTINE geo_coord_read_unit
423
424
425!> Legge da un'unità di file il contenuto dell'oggetto \a this.
426!! Il record da leggere deve essere stato scritto con la ::write_unit
427!! e, nel caso \a this sia un vettore, la lunghezza del record e quella
428!! del vettore devono essere accordate. Il metodo controlla se il file è
429!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
430SUBROUTINE geo_coord_vect_read_unit(this, unit)
431TYPE(geo_coord) :: this(:) !< oggetto da leggere
432INTEGER, INTENT(in) :: unit !< unità da cui leggere
433
434CHARACTER(len=40) :: form
435INTEGER :: i
436
437INQUIRE(unit, form=form)
438IF (form == 'FORMATTED') THEN
439 read(unit,*) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
440!TODO bug gfortran compiler !
441!missing values are unredeable when formatted
442ELSE
443 READ(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
444ENDIF
445
446END SUBROUTINE geo_coord_vect_read_unit
447
448
449!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
450!! Il record scritto potrà successivamente essere letto con la ::read_unit.
451!! Il metodo controlla se il file è
452!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
453SUBROUTINE geo_coord_write_unit(this, unit)
454TYPE(geo_coord),INTENT(in) :: this !< oggetto da scrivere
455INTEGER, INTENT(in) :: unit !< unità su cui scrivere
456
457CALL geo_coord_vect_write_unit((/this/), unit)
458
459END SUBROUTINE geo_coord_write_unit
461
462!> Scrive su un'unità di file il contenuto dell'oggetto \a this.
463!! Il record scritto potrà successivamente essere letto con la ::read_unit.
464!! Il metodo controlla se il file è
465!! aperto per un I/O formattato o non formattato e fa la cosa giusta.
466SUBROUTINE geo_coord_vect_write_unit(this, unit)
467TYPE(geo_coord),INTENT(in) :: this(:) !< oggetto da scrivere
468INTEGER, INTENT(in) :: unit !< unità su cui scrivere
469
470CHARACTER(len=40) :: form
471INTEGER :: i
472
473INQUIRE(unit, form=form)
474IF (form == 'FORMATTED') THEN
475 WRITE(unit,*) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
476!TODO bug gfortran compiler !
477!missing values are unredeable when formatted
478ELSE
479 WRITE(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
480ENDIF
481
482END SUBROUTINE geo_coord_vect_write_unit
483
484
485!> Restituisce la distanza in m tra 2 oggetti geo_coord.
486!! La distanza è calcolata approssimativamente ed è valida per piccoli angoli.
487FUNCTION geo_coord_dist(this, that) RESULT(dist)
489TYPE(geo_coord), INTENT (IN) :: this !< primo punto
490TYPE(geo_coord), INTENT (IN) :: that !< secondo punto
491REAL(kind=fp_geo) :: dist !< distanza in metri
492
493REAL(kind=fp_geo) :: x,y
494! Distanza approssimata, valida per piccoli angoli
495
496x = (getlon(this)-getlon(that))*cos(((getlat(this)+getlat(that))/2.)*degrad)
497y = getlat(this)-getlat(that)
498dist = sqrt(x**2 + y**2)*degrad*rearth
499
500END FUNCTION geo_coord_dist
501
502
503!> Determina se il punto indicato da \a this è contenuto in un rettangolo.
504!! Il rettangolo è orientato parallelamente agli assi del sistema,
505!! i suoi vertici sud-ovest e nord-est sono specificati da altri due punti.
506!! La funzione restituisce \c .TRUE. anche se il punto si trova sulla frontiera
507!! del rettangolo.
508!! Tutti gli oggetti devono essere già stati
509!! convertiti ad un sistema di coordinate comune, altrimenti viene
510!! restituito \c .FALSE. .
511FUNCTION geo_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
512TYPE(geo_coord),INTENT(IN) :: this !< oggetto di cui determinare la posizione
513TYPE(geo_coord),INTENT(IN) :: coordmin !< vertice sud-ovest del rettangolo
514TYPE(geo_coord),INTENT(IN) :: coordmax !< vertice nord-est del rettangolo
515LOGICAL :: res !< \c .TRUE. se \a this è dentro il rettangolo o sul bordo e \c .FALSE. se è fuori
516
517res = (geo_coord_ure(this, coordmin) .AND. geo_coord_lle(this,coordmax))
518
519END FUNCTION geo_coord_inside_rectang
520
521
522! ===================
523! == geo_coordvect ==
524! ===================
525!> Costruisce un oggetto \a geo_coordvect con i parametri opzionali forniti.
526!! Se sono presenti \a lon e \a lat, inizializza le coordinate geografiche
527!! ignorando \a utme e \a utmn, mentre se sono specificati \a utme e \a utmn
528!! succede il contrario; non è possibile specificare le coordinate in entrambi
529!! i sistemi, usare eventualmente \a to_geo. Se non viene passato nessun parametro
530!! opzionale l'oggetto è inizializzato a valore mancante.
531!! Il numero di punti dell'oggetto finale sarà uguale all'estensione
532!! del più breve vettore della coppia fornita.
533RECURSIVE SUBROUTINE geo_coordvect_init(this, lon, lat)
534TYPE(geo_coordvect), INTENT(OUT) :: this !< oggetto da inizializzare
535REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon(:) !< longitudine geografica
536REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat(:) !< latitudine geografica
537
538IF (PRESENT(lon) .AND. PRESENT(lat)) THEN
539 this%vsize = min(SIZE(lon), SIZE(lat))
540 ALLOCATE(this%ll(this%vsize,2))
541 this%ll(1:this%vsize,1) = lon(1:this%vsize)
542 this%ll(1:this%vsize,2) = lat(1:this%vsize)
543ELSE
544 this%vsize = 0
545 NULLIFY(this%ll)
546ENDIF
547this%vtype = 0 !?
549END SUBROUTINE geo_coordvect_init
550
551
552!> Distrugge l'oggetto in maniera pulita, liberando l'eventuale spazio
553!! dinamicamente allocato.
554SUBROUTINE geo_coordvect_delete(this)
555TYPE(geo_coordvect), INTENT(INOUT) :: this
556
557IF (ASSOCIATED(this%ll)) DEALLOCATE(this%ll)
558this%vsize = 0
559this%vtype = 0
560
561END SUBROUTINE geo_coordvect_delete
562
563
564!> Restituisce il valore di uno o più componenti di un oggetto \a geo_coordvect.
565!! Qualsiasi combinazione dei parametri opzionali è consentita; se
566!! il tipo di coordinata richiesta non è stato inizializzato né calcolato,
567!! restituisce il corrispondente valore mancante.
568!! Se forniti, i parametri \a lon, \a lat, \a utme, \a utmn devono essere
569!! dichiarati come puntatori che vengono
570!! allocati dalla \a getval stessa e che devono poi essere
571!! deallocati esplicitamente dal programma chiamante.
572SUBROUTINE geo_coordvect_getval(this, lon, lat)
573TYPE(geo_coordvect),INTENT(IN) :: this !< oggetto di cui restituire i componenti
574REAL(kind=fp_geo), OPTIONAL, POINTER :: lon(:) !< longitudine geografica
575REAL(kind=fp_geo), OPTIONAL, POINTER :: lat(:) !< latitudine geografica
576
577IF (PRESENT(lon)) THEN
578 IF (ASSOCIATED(this%ll)) THEN
579 ALLOCATE(lon(this%vsize))
580 lon(:) = this%ll(1:this%vsize,1)
581 ENDIF
582ENDIF
583IF (PRESENT(lat)) THEN
584 IF (ASSOCIATED(this%ll)) THEN
585 ALLOCATE(lat(this%vsize))
586 lat(:) = this%ll(1:this%vsize,2)
587 ENDIF
588ENDIF
589
590END SUBROUTINE geo_coordvect_getval
591
592
593SUBROUTINE geo_coordvect_import(this, unitsim, &
594#ifdef HAVE_SHAPELIB
595 shphandle, &
596#endif
597 nshp)
598TYPE(geo_coordvect), INTENT(OUT) :: this
599INTEGER,OPTIONAL,INTENT(IN) :: unitsim
600#ifdef HAVE_SHAPELIB
601TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
602#endif
603INTEGER,OPTIONAL,INTENT(IN) :: nshp
605REAL(kind=fp_geo),ALLOCATABLE :: llon(:), llat(:)
606REAL(kind=fp_geo) :: lv1,lv2,lv3,lv4
607INTEGER :: i, lvsize
608CHARACTER(len=40) :: lname
609#ifdef HAVE_SHAPELIB
610TYPE(shpobject) :: shpobj
611#endif
612
613IF (PRESENT(unitsim)) THEN
614 ! Leggo l'intestazione
615 READ(unitsim,*,END=10)lvsize,lv1,lv2,lv3,lv4,lname
616 ALLOCATE(llon(lvsize+1), llat(lvsize+1))
617 ! Leggo il poligono
618 READ(unitsim,*)(llon(i),llat(i), i=1,lvsize)
619 ! Lo chiudo se necessario
620 IF (llon(1) /= llon(lvsize) .OR. llat(1) /= llat(lvsize)) THEN
621 lvsize = lvsize + 1
622 llon(lvsize) = llon(1)
623 llat(lvsize) = llat(1)
624 ENDIF
625 ! Lo inserisco nel mio oggetto
626 CALL init(this, lon=llon(1:lvsize), lat=llat(1:lvsize))
627 this%vtype = geo_coordvect_polygon ! Sempre un poligono
628
629 DEALLOCATE(llon, llat)
630 RETURN
63110 CALL raise_error('nella lettura del file '//trim(to_char(unitsim)))
632 DEALLOCATE(llon, llat) ! End of file, ritorno un oggetto non assegnato
633#ifdef HAVE_SHAPELIB
634ELSE IF (PRESENT(shphandle) .AND. PRESENT(nshp)) THEN
635 ! Leggo l'oggetto shape
636 shpobj = shpreadobject(shphandle, nshp)
637 IF (.NOT.shpisnull(shpobj)) THEN
638 ! Lo inserisco nel mio oggetto
639 CALL init(this, lon=real(shpobj%padfx,kind=fp_geo), &
640 lat=real(shpobj%padfy,kind=fp_geo))
641 this%vtype = shpobj%nshptype
642 CALL shpdestroyobject(shpobj)
643 ELSE
644 CALL init(this)
645 ENDIF
646#endif
647ENDIF
648
649END SUBROUTINE geo_coordvect_import
650
651
652SUBROUTINE geo_coordvect_export(this, unitsim, &
653#ifdef HAVE_SHAPELIB
654 shphandle, &
655#endif
656 nshp)
657TYPE(geo_coordvect), INTENT(INOUT) :: this
658INTEGER,OPTIONAL,INTENT(IN) :: unitsim
659#ifdef HAVE_SHAPELIB
660TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
661#endif
662INTEGER,OPTIONAL,INTENT(IN) :: nshp
663
664INTEGER :: i, lnshp
665#ifdef HAVE_SHAPELIB
666TYPE(shpobject) :: shpobj
667#endif
668
669IF (PRESENT(unitsim)) THEN
670 IF (this%vsize > 0) THEN
671 ! Scrivo l'intestazione
672 WRITE(unitsim,*)SIZE(this%ll,1),-1.,5000.,-0.1,1.1,'Area'
673 ! Scrivo il poligono
674 WRITE(unitsim,*)(this%ll(i,1:2), i=1,this%vsize)
675 ELSE
676 CALL raise_warning('oggetto geo_coordvect vuoto, non scrivo niente in '// &
677 trim(to_char(unitsim)))
678 ENDIF
679#ifdef HAVE_SHAPELIB
680ELSE IF (PRESENT(shphandle)) THEN
681 IF (PRESENT(nshp)) THEN
682 lnshp = nshp
683 ELSE
684 lnshp = -1 ! -1 = append
685 ENDIF
686 ! Creo l'oggetto shape inizializzandolo con il mio oggetto
687 shpobj = shpcreatesimpleobject(this%vtype, this%vsize, &
688 REAL(this%ll(1:this%vsize,1),kind=fp_d), &
689 REAL(this%ll(1:this%vsize,2),kind=fp_d))
690 IF (.NOT.shpisnull(shpobj)) THEN
691 ! Lo scrivo nello shapefile
692 i=shpwriteobject(shphandle, lnshp, shpobj)
693 CALL shpdestroyobject(shpobj)
694 ENDIF
695#endif
696ENDIF
697
698END SUBROUTINE geo_coordvect_export
700!> Importa un vettore di oggetti \a geo_coordvect da un file in
701!! formato testo o in formato \a shapefile.
702!! Il parametro \a this è un puntatore che sarà allocato a cura del metodo stesso e
703!! dovrà invece essere deallocato da parte del programma chiamante
704!! dopo aver chiamato il metodo \a delete per ogni suo elemento.
705!! In caso di errore nella fase iniziale di importazione,
706!! \a this non verrà associato, e quindi è opportuno testare
707!! \code
708!! IF (ASSOCIATED(my_coord_vect)) THEN...
709!! \endcode
710!! nel programma chiamante
711!! per intrappolare eventuale condizioni di errore (tipicamente file non
712!! trovato o in un formato non compatibile).
713!! Entrambi i formati di ingresso non contengono informazioni sul tipo
714!! di coordinate dei dati
715!! (per il formato shapefile è possibile solo con delle estensioni non standard),
716!! per cui questa informazione, se desiderata, deve essere fornita dal
717!! programma chiamante.
718SUBROUTINE geo_coordvect_importvect(this, shpfilesim, shpfile)
719TYPE(geo_coordvect),POINTER :: this(:) !< puntatore all'oggetto su cui importare i dati, viene allocato dalla \a import stessa
720CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim !< nome del file in formato testo "SIM", il parametro deve essere fornito solo se si vuole importare da un file di quel tipo
721CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile !< nome delllo shapefile, il parametro deve essere fornito solo se si vuole importare da un file di quel tipo
722
723REAL(kind=fp_geo) :: inu
724REAL(kind=fp_d) :: minb(4), maxb(4)
725INTEGER :: i, u, ns, lvsize, shptype, dbfnf, dbfnr
726CHARACTER(len=40) :: lname
727#ifdef HAVE_SHAPELIB
728TYPE(shpfileobject) :: shphandle
729#endif
730
731NULLIFY(this)
732
733IF (PRESENT(shpfilesim)) THEN
734 u = getunit()
735 OPEN(u, file=shpfilesim, status='old', err=30)
736 ns = 0 ! Conto il numero di shape contenute
737 DO WHILE(.true.)
738 READ(u,*,END=10,ERR=20)lvsize,inu,inu,inu,inu,lname
739 READ(u,*,END=20,ERR=20)(inu,inu, i=1,lvsize)
740 ns = ns + 1
741 ENDDO
74210 CONTINUE
743 IF (ns > 0) THEN ! Alloco e leggo il mio oggetto
744 ALLOCATE(this(ns))
745 rewind(u)
746 DO i = 1, ns
747 CALL import(this(i), unitsim=u)
748 ENDDO
749 ENDIF
75020 CONTINUE
751 CLOSE(u)
752 IF (.NOT.ASSOCIATED(this)) THEN
753 CALL raise_warning('file '//trim(shpfilesim)//' vuoto o corrotto')
754 ENDIF
755 RETURN
75630 CONTINUE
757 CALL raise_error('Impossibile aprire il file '//trim(shpfile))
758 RETURN
759
760ELSE IF (PRESENT(shpfile)) THEN
761#ifdef HAVE_SHAPELIB
762 shphandle = shpopen(trim(shpfile), 'rb')
763 IF (shpfileisnull(shphandle)) THEN
764 CALL raise_error('Impossibile aprire lo shapefile '//trim(shpfile))
765 RETURN
766 ENDIF
767 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr) ! Ottengo le info sul file
768 IF (ns > 0) THEN ! Alloco e leggo il mio oggetto
769 ALLOCATE(this(ns))
770 this(:)%vtype = shptype
771 DO i = 1, ns
772 CALL import(this(i), shphandle=shphandle, nshp=i-1)
773 ENDDO
774 ENDIF
775 CALL shpclose(shphandle)
776 RETURN
777#endif
778ENDIF
779
780END SUBROUTINE geo_coordvect_importvect
781
782
783!> Esporta un vettore di oggetti \a geo_coordvect su un file in
784!! formato testo o in formato \a shapefile.
785SUBROUTINE geo_coordvect_exportvect(this, shpfilesim, shpfile, append)
786TYPE(geo_coordvect) :: this(:) !< oggetto da esportare
787CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim !< nome del file in formato testo "SIM", il parametro deve essere fornito solo se si vuole esportare su un file di quel tipo
788CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile !< nome dello shapefile, il parametro deve essere fornito solo se si vuole esportare su un file di quel tipo
789!> sistema di coordinate (proiezione) dei dati,
790!! usare i parametri \a ::proj_geo (default) o \a ::proj_utm,
791!! ha senso se \a this, a seguito di una chiamata a \a ::to_geo o a \a ::to_utm,
792!! contiene le coordinate in entrambi i sistemi, altrimenti i dati vengono
793!! esportati automaticamente nel solo sistema disponibile
794LOGICAL,INTENT(in),OPTIONAL :: append !< se è presente e vale \c .TRUE. , ::export accoda all'eventuale file esistente anziché sovrascriverlo
795
796REAL(kind=fp_d) :: minb(4), maxb(4)
797INTEGER :: i, u, ns, shptype, dbfnf, dbfnr
798LOGICAL :: lappend
799#ifdef HAVE_SHAPELIB
800TYPE(shpfileobject) :: shphandle
801#endif
802
803IF (PRESENT(append)) THEN
804 lappend = append
805ELSE
806 lappend = .false.
807ENDIF
808IF (PRESENT(shpfilesim)) THEN
809 u = getunit()
810 IF (lappend) THEN
811 OPEN(u, file=shpfilesim, status='unknown', position='append', err=30)
812 ELSE
813 OPEN(u, file=shpfilesim, status='unknown', err=30)
814 ENDIF
815 DO i = 1, SIZE(this)
816 CALL export(this(i), unitsim=u)
817 ENDDO
818 CLOSE(u)
819 RETURN
82030 CONTINUE
821 CALL raise_error('Impossibile aprire il file '//trim(shpfile))
822 RETURN
823ELSE IF (PRESENT(shpfile)) THEN
824#ifdef HAVE_SHAPELIB
825 IF (lappend) THEN
826 shphandle = shpopen(trim(shpfile), 'r+b')
827 IF (shpfileisnull(shphandle)) THEN ! shpopen funziona solo su file esistenti
828 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
829 ENDIF
830 ELSE
831 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
832 ENDIF
833 IF (shpfileisnull(shphandle)) THEN
834 CALL raise_error('Impossibile aprire lo shapefile '//trim(shpfile))
835 RETURN
836 ENDIF
837 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr) ! Ottengo le info sul file
838 DO i = 1, SIZE(this)
839 IF (i > ns .OR. lappend) THEN ! Append shape
840 CALL export(this(i), shphandle=shphandle)
841 ELSE ! Overwrite shape
842 CALL export(this(i), shphandle=shphandle, nshp=i-1)
843 ENDIF
844 ENDDO
845 CALL shpclose(shphandle)
846 RETURN
847#endif
848ENDIF
849
850END SUBROUTINE geo_coordvect_exportvect
851
852
853!> Determina se il punto indicato da \a this si trova
854!! dentro o fuori dal poligono descritto dall'oggetto \a poly.
855!! Funziona anche se la topologia di \a poly non è poligonale,
856!! forzandone la chiusura; usa un algoritmo di ricerca del numero di
857!! intersezioni, come indicato in
858!! comp.graphics.algorithms FAQ (http://www.faqs.org/faqs/graphics/algorithms-faq/)
859!! o in
860!! http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
861FUNCTION geo_coord_inside(this, poly) RESULT(inside)
862TYPE(geo_coord), INTENT(IN) :: this !< oggetto di cui determinare la posizione
863TYPE(geo_coordvect), INTENT(IN) :: poly !< poligono contenitore
864LOGICAL :: inside
865
866INTEGER :: i, j, starti
867
868inside = .false.
869IF (all(poly%ll(1,:) == poly%ll(poly%vsize,:))) THEN ! Poligono chiuso
870 starti = 2
871 j = 1
872ELSE ! Poligono non chiuso
873 starti = 1
874 j = poly%vsize
875ENDIF
876DO i = starti, poly%vsize
877 IF ((poly%ll(i,2) <= getlat(this) .AND. &
878 getlat(this) < poly%ll(j,2)) .OR. &
879 (poly%ll(j,2) <= getlat(this) .AND. &
880 getlat(this) < poly%ll(i,2))) THEN
881 IF (getlon(this) < (poly%ll(j,1) - poly%ll(i,1)) * &
882 (getlat(this) - poly%ll(i,2)) / &
883 (poly%ll(j,2) - poly%ll(i,2)) + poly%ll(i,1)) THEN
884 inside = .NOT. inside
885 ENDIF
886 ENDIF
887 j = i
888ENDDO
889
890END FUNCTION geo_coord_inside
891
892
893ELEMENTAL FUNCTION c_e_geo_coord(this) result (res)
894TYPE(geo_coord),INTENT(in) :: this
895LOGICAL :: res
896
897res = .not. this == geo_coord_miss
898
899end FUNCTION c_e_geo_coord
900
901
902character(len=80) function to_char_geo_coord(this)
903TYPE(geo_coord),INTENT(in) :: this
904
905to_char_geo_coord = "GEO_COORD: Lon="// &
906 trim(to_char(getlon(this),miss="Missing lon",form="(f11.5)"))//&
907 " Lat="// &
908 trim(to_char(getlat(this),miss="Missing lat",form="(f11.5)"))
909
910end function to_char_geo_coord
911
912
913subroutine display_geo_coord(this)
914TYPE(geo_coord),INTENT(in) :: this
915
916print*,trim(to_char(this))
917
918end subroutine display_geo_coord
919
920
921END MODULE geo_coord_class
Detructors for the two classes.
Export one or more geo_coordvect objects to a plain text file or to a file in ESRI/Shapefile format.
Methods for returning the value of object members.
Import one or more geo_coordvect objects from a plain text file or for a file in ESRI/Shapefile forma...
Constructors for the two classes.
Determine whether a point lies inside a polygon or a rectangle.
Read a single geo_coord object or an array of geo_coord objects from a Fortran FORMATTED or UNFORMATT...
Represent geo_coord object in a pretty string.
Write a single geo_coord object or an array of geo_coord objects to a Fortran FORMATTED or UNFORMATTE...
Utilities for CHARACTER variables.
Costanti fisiche (DOUBLEPRECISION).
Gestione degli errori.
Utilities for managing files.
Classes for handling georeferenced sparse points in geographical corodinates.
Definition of constants to be used for declaring variables of a desired type.
Definition kinds.F90:245
Definitions of constants and functions for working with missing values.
Derived type defining an isolated georeferenced point on Earth in polar geographical coordinates.
Derived type defining a one-dimensional array of georeferenced points with an associated topology (is...

Generated with Doxygen.