libsim Versione 7.2.4
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
249
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
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
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
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
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
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
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
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
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
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
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 !?
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
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
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
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.