libsim Versione 7.1.11
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
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
50INTEGER, PARAMETER :: fp_geo=fp_d
51real(kind=fp_geo), PARAMETER :: fp_geo_miss=dmiss
52
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
70 PRIVATE
71 REAL(kind=fp_geo),POINTER :: ll(:,:) => null()
72 INTEGER :: vsize, vtype
73END TYPE geo_coordvect
74
76TYPE(geo_coord),PARAMETER :: geo_coord_miss=geo_coord(imiss,imiss)
77
78INTEGER, PARAMETER :: geo_coordvect_point = 1
79INTEGER, PARAMETER :: geo_coordvect_arc = 3
80INTEGER, PARAMETER :: geo_coordvect_polygon = 5
81INTEGER, PARAMETER :: geo_coordvect_multipoint = 8
82
83
87INTERFACE init
88 MODULE PROCEDURE geo_coord_init, geo_coordvect_init
89END INTERFACE
90
94INTERFACE delete
95 MODULE PROCEDURE geo_coord_delete, geo_coordvect_delete
96END INTERFACE
97
99INTERFACE getval
100 MODULE PROCEDURE geo_coord_getval, geo_coordvect_getval
101END INTERFACE
102
104INTERFACE OPERATOR (==)
105 MODULE PROCEDURE geo_coord_eq
106END INTERFACE
107
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
131INTERFACE import
132 MODULE PROCEDURE geo_coordvect_import, geo_coordvect_importvect
133END INTERFACE
134
137INTERFACE export
138 MODULE PROCEDURE geo_coordvect_export, geo_coordvect_exportvect
139END INTERFACE
140
143INTERFACE read_unit
144 MODULE PROCEDURE geo_coord_read_unit, geo_coord_vect_read_unit
145END INTERFACE
146
149INTERFACE write_unit
150 MODULE PROCEDURE geo_coord_write_unit, geo_coord_vect_write_unit
151END INTERFACE
152
154INTERFACE inside
155 MODULE PROCEDURE geo_coord_inside, geo_coord_inside_rectang
156END INTERFACE
157
159INTERFACE c_e
160 MODULE PROCEDURE c_e_geo_coord
161END INTERFACE
162
164INTERFACE to_char
165 MODULE PROCEDURE to_char_geo_coord
166END INTERFACE
167
169INTERFACE display
170 MODULE PROCEDURE display_geo_coord
171END INTERFACE
172
173CONTAINS
174
175
176! ===================
177! == geo_coord ==
178! ===================
185SUBROUTINE geo_coord_init(this, lon, lat, ilon, ilat)
186TYPE(geo_coord) :: this
187REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon
188REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat
189integer(kind=int_l), INTENT(IN), OPTIONAL :: ilon
190integer(kind=int_l), INTENT(IN), OPTIONAL :: ilat
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
210SUBROUTINE geo_coord_delete(this)
211TYPE(geo_coord), INTENT(INOUT) :: this
212
213this%ilon = imiss
214this%ilat = imiss
215
216END SUBROUTINE geo_coord_delete
217
222elemental SUBROUTINE geo_coord_getval(this, lon, lat,ilon,ilat)
223TYPE(geo_coord),INTENT(IN) :: this
224REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lon
225REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lat
226integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilon
227integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilat
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
242elemental FUNCTION getilat(this)
243TYPE(geo_coord),INTENT(IN) :: this
244integer(kind=int_l) :: getilat
245
246getilat = this%ilat
247
248END FUNCTION getilat
254elemental FUNCTION getlat(this)
255TYPE(geo_coord),INTENT(IN) :: this
256real(kind=fp_geo) :: getlat
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
267
272elemental FUNCTION getilon(this)
273TYPE(geo_coord),INTENT(IN) :: this
274integer(kind=int_l) :: getilon
276getilon = this%ilon
277
278END FUNCTION getilon
279
280
285elemental FUNCTION getlon(this)
286TYPE(geo_coord),INTENT(IN) :: this
287real(kind=fp_geo) :: getlon
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
299
300elemental FUNCTION geo_coord_eq(this, that) RESULT(res)
301TYPE(geo_coord),INTENT(IN) :: this, that
302LOGICAL :: res
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
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
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
344elemental FUNCTION geo_coord_lt(this, that) RESULT(res)
345TYPE(geo_coord),INTENT(IN) :: this, that
346LOGICAL :: res
347
348res = this%ilon < that%ilon
349
350if ( this%ilon == that%ilon ) then
351 res = this%ilat < that%ilat
352end if
354
355END FUNCTION geo_coord_lt
356
357
360elemental FUNCTION geo_coord_le(this, that) RESULT(res)
361TYPE(geo_coord),INTENT(IN) :: this, that
362LOGICAL :: res
364res = .not. this > that
365
366END FUNCTION geo_coord_le
367
368
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
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
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
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
410
416SUBROUTINE geo_coord_read_unit(this, unit)
417TYPE(geo_coord),INTENT(out) :: this
418INTEGER, INTENT(in) :: unit
419
420CALL geo_coord_vect_read_unit((/this/), unit)
421
422END SUBROUTINE geo_coord_read_unit
423
424
430SUBROUTINE geo_coord_vect_read_unit(this, unit)
431TYPE(geo_coord) :: this(:)
432INTEGER, INTENT(in) :: unit
433
434CHARACTER(len=40) :: form
435INTEGER :: i
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
453SUBROUTINE geo_coord_write_unit(this, unit)
454TYPE(geo_coord),INTENT(in) :: this
455INTEGER, INTENT(in) :: unit
456
457CALL geo_coord_vect_write_unit((/this/), unit)
458
459END SUBROUTINE geo_coord_write_unit
460
461
466SUBROUTINE geo_coord_vect_write_unit(this, unit)
467TYPE(geo_coord),INTENT(in) :: this(:)
468INTEGER, INTENT(in) :: unit
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
487FUNCTION geo_coord_dist(this, that) RESULT(dist)
489TYPE(geo_coord), INTENT (IN) :: this
490TYPE(geo_coord), INTENT (IN) :: that
491REAL(kind=fp_geo) :: dist
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
511FUNCTION geo_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
512TYPE(geo_coord),INTENT(IN) :: this
513TYPE(geo_coord),INTENT(IN) :: coordmin
514TYPE(geo_coord),INTENT(IN) :: coordmax
515LOGICAL :: res
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! ===================
533RECURSIVE SUBROUTINE geo_coordvect_init(this, lon, lat)
534TYPE(geo_coordvect), INTENT(OUT) :: this
535REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon(:)
536REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat(:)
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 !?
548
549END SUBROUTINE geo_coordvect_init
550
551
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
572SUBROUTINE geo_coordvect_getval(this, lon, lat)
573TYPE(geo_coordvect),INTENT(IN) :: this
574REAL(kind=fp_geo), OPTIONAL, POINTER :: lon(:)
575REAL(kind=fp_geo), OPTIONAL, POINTER :: lat(:)
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
604
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
699
718SUBROUTINE geo_coordvect_importvect(this, shpfilesim, shpfile)
719TYPE(geo_coordvect),POINTER :: this(:)
720CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim
721CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile
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
785SUBROUTINE geo_coordvect_exportvect(this, shpfilesim, shpfile, append)
786TYPE(geo_coordvect) :: this(:)
787CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim
788CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile
794LOGICAL,INTENT(in),OPTIONAL :: append
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
861FUNCTION geo_coord_inside(this, poly) RESULT(inside)
862TYPE(geo_coord), INTENT(IN) :: this
863TYPE(geo_coordvect), INTENT(IN) :: poly
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
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).
Definition: phys_const.f90:72
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:251
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.