libsim  Versione7.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 
28 MODULE geo_coord_class
29 USE kinds
30 USE err_handling
33 !USE doubleprecision_phys_const
35 #ifdef HAVE_SHAPELIB
36 USE shapelib
37 !, ONLY: shpobject, shpreadobject, shpdestroyobject, shpwriteobject, &
38 ! shpcreatesimpleobject
39 #endif
40 IMPLICIT NONE
41 
42 
50 INTEGER, PARAMETER :: fp_geo=fp_d
51 real(kind=fp_geo), PARAMETER :: fp_geo_miss=dmiss
52 
55 TYPE geo_coord
56  PRIVATE
57  INTEGER(kind=int_l) :: ilon, ilat
58 END TYPE geo_coord
59 
60 TYPE geo_coord_r
61  PRIVATE
62  REAL(kind=fp_geo) :: lon, lat
63 END TYPE geo_coord_r
64 
65 
69 TYPE geo_coordvect
70  PRIVATE
71  REAL(kind=fp_geo),POINTER :: ll(:,:) => null()
72  INTEGER :: vsize, vtype
73 END TYPE geo_coordvect
74 
76 TYPE(geo_coord),PARAMETER :: geo_coord_miss=geo_coord(imiss,imiss)
77 
78 INTEGER, PARAMETER :: geo_coordvect_point = 1
79 INTEGER, PARAMETER :: geo_coordvect_arc = 3
80 INTEGER, PARAMETER :: geo_coordvect_polygon = 5
81 INTEGER, PARAMETER :: geo_coordvect_multipoint = 8
82 
83 
87 INTERFACE init
88  MODULE PROCEDURE geo_coord_init, geo_coordvect_init
89 END INTERFACE
90 
94 INTERFACE delete
95  MODULE PROCEDURE geo_coord_delete, geo_coordvect_delete
96 END INTERFACE
97 
99 INTERFACE getval
100  MODULE PROCEDURE geo_coord_getval, geo_coordvect_getval
101 END INTERFACE
102 
104 INTERFACE OPERATOR (==)
105  MODULE PROCEDURE geo_coord_eq
106 END INTERFACE
107 
109 INTERFACE OPERATOR (/=)
110  MODULE PROCEDURE geo_coord_ne
111 END INTERFACE
112 
113 INTERFACE OPERATOR (>=)
114  MODULE PROCEDURE geo_coord_ge
115 END INTERFACE
116 
117 INTERFACE OPERATOR (>)
118  MODULE PROCEDURE geo_coord_gt
119 END INTERFACE
120 
121 INTERFACE OPERATOR (<=)
122  MODULE PROCEDURE geo_coord_le
123 END INTERFACE
124 
125 INTERFACE OPERATOR (<)
126  MODULE PROCEDURE geo_coord_lt
127 END INTERFACE
128 
131 INTERFACE import
132  MODULE PROCEDURE geo_coordvect_import, geo_coordvect_importvect
133 END INTERFACE
134 
137 INTERFACE export
138  MODULE PROCEDURE geo_coordvect_export, geo_coordvect_exportvect
139 END INTERFACE
140 
143 INTERFACE read_unit
144  MODULE PROCEDURE geo_coord_read_unit, geo_coord_vect_read_unit
145 END INTERFACE
146 
149 INTERFACE write_unit
150  MODULE PROCEDURE geo_coord_write_unit, geo_coord_vect_write_unit
151 END INTERFACE
152 
154 INTERFACE inside
155  MODULE PROCEDURE geo_coord_inside, geo_coord_inside_rectang
156 END INTERFACE
157 
159 INTERFACE c_e
160  MODULE PROCEDURE c_e_geo_coord
161 END INTERFACE
162 
164 INTERFACE to_char
165  MODULE PROCEDURE to_char_geo_coord
166 END INTERFACE
167 
169 INTERFACE display
170  MODULE PROCEDURE display_geo_coord
171 END INTERFACE
172 
173 CONTAINS
174 
175 
176 ! ===================
177 ! == geo_coord ==
178 ! ===================
185 SUBROUTINE geo_coord_init(this, lon, lat, ilon, ilat)
186 TYPE(geo_coord) :: this
187 REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon
188 REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat
189 integer(kind=int_l), INTENT(IN), OPTIONAL :: ilon
190 integer(kind=int_l), INTENT(IN), OPTIONAL :: ilat
191 
192 real(kind=fp_geo) :: llon,llat
193 
194 CALL optio(ilon, this%ilon)
195 CALL optio(ilat, this%ilat)
196 
197 if (.not. c_e(this%ilon)) then
198  CALL optio(lon, llon)
199  if (c_e(llon)) this%ilon=nint(llon*1.d5)
200 end if
201 
202 if (.not. c_e(this%ilat)) then
203  CALL optio(lat, llat)
204  if (c_e(llat)) this%ilat=nint(llat*1.d5)
205 end if
206 
207 END SUBROUTINE geo_coord_init
208 
210 SUBROUTINE geo_coord_delete(this)
211 TYPE(geo_coord), INTENT(INOUT) :: this
212 
213 this%ilon = imiss
214 this%ilat = imiss
215 
216 END SUBROUTINE geo_coord_delete
217 
222 elemental SUBROUTINE geo_coord_getval(this, lon, lat,ilon,ilat)
223 TYPE(geo_coord),INTENT(IN) :: this
224 REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lon
225 REAL(kind=fp_geo), INTENT(OUT), OPTIONAL :: lat
226 integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilon
227 integer(kind=int_l), INTENT(OUT), OPTIONAL :: ilat
228 
229 IF (PRESENT(ilon)) ilon = getilon(this)
230 IF (PRESENT(ilat)) ilat = getilat(this)
231 
232 IF (PRESENT(lon)) lon = getlon(this)
233 IF (PRESENT(lat)) lat = getlat(this)
234 
235 END SUBROUTINE geo_coord_getval
236 
237 
242 elemental FUNCTION getilat(this)
243 TYPE(geo_coord),INTENT(IN) :: this
244 integer(kind=int_l) :: getilat
245 
246 getilat = this%ilat
247 
248 END FUNCTION getilat
249 
254 elemental FUNCTION getlat(this)
255 TYPE(geo_coord),INTENT(IN) :: this
256 real(kind=fp_geo) :: getlat
257 integer(kind=int_l) :: ilat
258 
259 ilat=getilat(this)
260 if (c_e(ilat)) then
261  getlat = ilat*1.d-5
262 else
263  getlat=fp_geo_miss
264 end if
265 
266 END FUNCTION getlat
272 elemental FUNCTION getilon(this)
273 TYPE(geo_coord),INTENT(IN) :: this
274 integer(kind=int_l) :: getilon
275 
276 getilon = this%ilon
278 END FUNCTION getilon
280 
285 elemental FUNCTION getlon(this)
286 TYPE(geo_coord),INTENT(IN) :: this
287 real(kind=fp_geo) :: getlon
288 integer(kind=int_l) :: ilon
289 
290 ilon=getilon(this)
291 if (c_e(ilon)) then
292  getlon = ilon*1.d-5
293 else
294  getlon=fp_geo_miss
295 end if
296 
297 END FUNCTION getlon
298 
299 
300 elemental FUNCTION geo_coord_eq(this, that) RESULT(res)
301 TYPE(geo_coord),INTENT(IN) :: this, that
302 LOGICAL :: res
303 
304 res = (this%ilon == that%ilon .AND. this%ilat == that%ilat)
305 
306 END FUNCTION geo_coord_eq
308 
309 elemental FUNCTION geo_coord_ne(this, that) RESULT(res)
310 TYPE(geo_coord),INTENT(IN) :: this, that
311 LOGICAL :: res
312 
313 res = .not. this == that
314 
315 END FUNCTION geo_coord_ne
316 
319 elemental FUNCTION geo_coord_gt(this, that) RESULT(res)
320 TYPE(geo_coord),INTENT(IN) :: this, that
321 LOGICAL :: res
322 
323 res = this%ilon > that%ilon
324 
325 if ( this%ilon == that%ilon ) then
326  res = this%ilat > that%ilat
327 end if
328 
329 END FUNCTION geo_coord_gt
330 
331 
334 elemental FUNCTION geo_coord_ge(this, that) RESULT(res)
335 TYPE(geo_coord),INTENT(IN) :: this, that
336 LOGICAL :: res
337 
338 res = .not. this < that
339 
340 END FUNCTION geo_coord_ge
344 elemental FUNCTION geo_coord_lt(this, that) RESULT(res)
345 TYPE(geo_coord),INTENT(IN) :: this, that
346 LOGICAL :: res
348 res = this%ilon < that%ilon
349 
350 if ( this%ilon == that%ilon ) then
351  res = this%ilat < that%ilat
352 end if
353 
354 
355 END FUNCTION geo_coord_lt
356 
360 elemental FUNCTION geo_coord_le(this, that) RESULT(res)
361 TYPE(geo_coord),INTENT(IN) :: this, that
362 LOGICAL :: res
363 
364 res = .not. this > that
365 
366 END FUNCTION geo_coord_le
368 
371 elemental FUNCTION geo_coord_ure(this, that) RESULT(res)
372 TYPE(geo_coord),INTENT(IN) :: this, that
373 LOGICAL :: res
374 
375 res = (this%ilon >= that%ilon .AND. this%ilat >= that%ilat)
376 
377 END FUNCTION geo_coord_ure
378 
381 elemental FUNCTION geo_coord_ur(this, that) RESULT(res)
382 TYPE(geo_coord),INTENT(IN) :: this, that
383 LOGICAL :: res
385 res = (this%ilon > that%ilon .AND. this%ilat > that%ilat)
386 
387 END FUNCTION geo_coord_ur
388 
389 
392 elemental FUNCTION geo_coord_lle(this, that) RESULT(res)
393 TYPE(geo_coord),INTENT(IN) :: this, that
394 LOGICAL :: res
395 
396 res = (this%ilon <= that%ilon .AND. this%ilat <= that%ilat)
397 
398 END FUNCTION geo_coord_lle
399 
402 elemental FUNCTION geo_coord_ll(this, that) RESULT(res)
403 TYPE(geo_coord),INTENT(IN) :: this, that
404 LOGICAL :: res
405 
406 res = (this%ilon < that%ilon .AND. this%ilat < that%ilat)
407 
408 END FUNCTION geo_coord_ll
410 
416 SUBROUTINE geo_coord_read_unit(this, unit)
417 TYPE(geo_coord),INTENT(out) :: this
418 INTEGER, INTENT(in) :: unit
419 
420 CALL geo_coord_vect_read_unit((/this/), unit)
422 END SUBROUTINE geo_coord_read_unit
423 
424 
430 SUBROUTINE geo_coord_vect_read_unit(this, unit)
431 TYPE(geo_coord) :: this(:)
432 INTEGER, INTENT(in) :: unit
433 
434 CHARACTER(len=40) :: form
435 INTEGER :: i
436 
437 INQUIRE(unit, form=form)
438 IF (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
442 ELSE
443  READ(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
444 ENDIF
445 
446 END SUBROUTINE geo_coord_vect_read_unit
447 
448 
453 SUBROUTINE geo_coord_write_unit(this, unit)
454 TYPE(geo_coord),INTENT(in) :: this
455 INTEGER, INTENT(in) :: unit
456 
457 CALL geo_coord_vect_write_unit((/this/), unit)
458 
459 END SUBROUTINE geo_coord_write_unit
460 
461 
466 SUBROUTINE geo_coord_vect_write_unit(this, unit)
467 TYPE(geo_coord),INTENT(in) :: this(:)
468 INTEGER, INTENT(in) :: unit
469 
470 CHARACTER(len=40) :: form
471 INTEGER :: i
472 
473 INQUIRE(unit, form=form)
474 IF (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
478 ELSE
479  WRITE(unit) (this(i)%ilon,this(i)%ilat, i=1,SIZE(this))
480 ENDIF
481 
482 END SUBROUTINE geo_coord_vect_write_unit
483 
487 FUNCTION geo_coord_dist(this, that) RESULT(dist)
489 TYPE(geo_coord), INTENT (IN) :: this
490 TYPE(geo_coord), INTENT (IN) :: that
491 REAL(kind=fp_geo) :: dist
492 
493 REAL(kind=fp_geo) :: x,y
494 ! Distanza approssimata, valida per piccoli angoli
495 
496 x = (getlon(this)-getlon(that))*cos(((getlat(this)+getlat(that))/2.)*degrad)
497 y = getlat(this)-getlat(that)
498 dist = sqrt(x**2 + y**2)*degrad*rearth
499 
500 END FUNCTION geo_coord_dist
501 
502 
511 FUNCTION geo_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
512 TYPE(geo_coord),INTENT(IN) :: this
513 TYPE(geo_coord),INTENT(IN) :: coordmin
514 TYPE(geo_coord),INTENT(IN) :: coordmax
515 LOGICAL :: res
516 
517 res = (geo_coord_ure(this, coordmin) .AND. geo_coord_lle(this,coordmax))
519 END FUNCTION geo_coord_inside_rectang
520 
521 
522 ! ===================
523 ! == geo_coordvect ==
524 ! ===================
533 RECURSIVE SUBROUTINE geo_coordvect_init(this, lon, lat)
534 TYPE(geo_coordvect), INTENT(OUT) :: this
535 REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lon(:)
536 REAL(kind=fp_geo), INTENT(IN), OPTIONAL :: lat(:)
537 
538 IF (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)
543 ELSE
544  this%vsize = 0
545  NULLIFY(this%ll)
546 ENDIF
547 this%vtype = 0 !?
548 
549 END SUBROUTINE geo_coordvect_init
550 
551 
554 SUBROUTINE geo_coordvect_delete(this)
555 TYPE(geo_coordvect), INTENT(INOUT) :: this
556 
557 IF (ASSOCIATED(this%ll)) DEALLOCATE(this%ll)
558 this%vsize = 0
559 this%vtype = 0
560 
561 END SUBROUTINE geo_coordvect_delete
562 
563 
572 SUBROUTINE geo_coordvect_getval(this, lon, lat)
573 TYPE(geo_coordvect),INTENT(IN) :: this
574 REAL(kind=fp_geo), OPTIONAL, POINTER :: lon(:)
575 REAL(kind=fp_geo), OPTIONAL, POINTER :: lat(:)
576 
577 IF (PRESENT(lon)) THEN
578  IF (ASSOCIATED(this%ll)) THEN
579  ALLOCATE(lon(this%vsize))
580  lon(:) = this%ll(1:this%vsize,1)
581  ENDIF
582 ENDIF
583 IF (PRESENT(lat)) THEN
584  IF (ASSOCIATED(this%ll)) THEN
585  ALLOCATE(lat(this%vsize))
586  lat(:) = this%ll(1:this%vsize,2)
587  ENDIF
588 ENDIF
589 
590 END SUBROUTINE geo_coordvect_getval
592 
593 SUBROUTINE geo_coordvect_import(this, unitsim, &
594 #ifdef HAVE_SHAPELIB
595  shphandle, &
596 #endif
597  nshp)
598 TYPE(geo_coordvect), INTENT(OUT) :: this
599 INTEGER,OPTIONAL,INTENT(IN) :: unitsim
600 #ifdef HAVE_SHAPELIB
601 TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
602 #endif
603 INTEGER,OPTIONAL,INTENT(IN) :: nshp
604 
605 REAL(kind=fp_geo),ALLOCATABLE :: llon(:), llat(:)
606 REAL(kind=fp_geo) :: lv1,lv2,lv3,lv4
607 INTEGER :: i, lvsize
608 CHARACTER(len=40) :: lname
609 #ifdef HAVE_SHAPELIB
610 TYPE(shpobject) :: shpobj
611 #endif
612 
613 IF (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
631 10 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
634 ELSE 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
647 ENDIF
648 
649 END SUBROUTINE geo_coordvect_import
650 
651 
652 SUBROUTINE geo_coordvect_export(this, unitsim, &
653 #ifdef HAVE_SHAPELIB
654  shphandle, &
655 #endif
656  nshp)
657 TYPE(geo_coordvect), INTENT(INOUT) :: this
658 INTEGER,OPTIONAL,INTENT(IN) :: unitsim
659 #ifdef HAVE_SHAPELIB
660 TYPE(shpfileobject),OPTIONAL,INTENT(INOUT) :: shphandle
661 #endif
662 INTEGER,OPTIONAL,INTENT(IN) :: nshp
663 
664 INTEGER :: i, lnshp
665 #ifdef HAVE_SHAPELIB
666 TYPE(shpobject) :: shpobj
667 #endif
668 
669 IF (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
680 ELSE 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
696 ENDIF
697 
698 END SUBROUTINE geo_coordvect_export
699 
718 SUBROUTINE geo_coordvect_importvect(this, shpfilesim, shpfile)
719 TYPE(geo_coordvect),POINTER :: this(:)
720 CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim
721 CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile
722 
723 REAL(kind=fp_geo) :: inu
724 REAL(kind=fp_d) :: minb(4), maxb(4)
725 INTEGER :: i, u, ns, lvsize, shptype, dbfnf, dbfnr
726 CHARACTER(len=40) :: lname
727 #ifdef HAVE_SHAPELIB
728 TYPE(shpfileobject) :: shphandle
729 #endif
730 
731 NULLIFY(this)
733 IF (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
742 10 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
750 20 CONTINUE
751  CLOSE(u)
752  IF (.NOT.ASSOCIATED(this)) THEN
753  CALL raise_warning('file '//trim(shpfilesim)//' vuoto o corrotto')
754  ENDIF
755  RETURN
756 30 CONTINUE
757  CALL raise_error('Impossibile aprire il file '//trim(shpfile))
758  RETURN
759 
760 ELSE 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
778 ENDIF
779 
780 END SUBROUTINE geo_coordvect_importvect
781 
782 
785 SUBROUTINE geo_coordvect_exportvect(this, shpfilesim, shpfile, append)
786 TYPE(geo_coordvect) :: this(:)
787 CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfilesim
788 CHARACTER(len=*),INTENT(in),OPTIONAL :: shpfile
794 LOGICAL,INTENT(in),OPTIONAL :: append
795 
796 REAL(kind=fp_d) :: minb(4), maxb(4)
797 INTEGER :: i, u, ns, shptype, dbfnf, dbfnr
798 LOGICAL :: lappend
799 #ifdef HAVE_SHAPELIB
800 TYPE(shpfileobject) :: shphandle
801 #endif
802 
803 IF (PRESENT(append)) THEN
804  lappend = append
805 ELSE
806  lappend = .false.
807 ENDIF
808 IF (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
820 30 CONTINUE
821  CALL raise_error('Impossibile aprire il file '//trim(shpfile))
822  RETURN
823 ELSE 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
848 ENDIF
849 
850 END SUBROUTINE geo_coordvect_exportvect
851 
852 
861 FUNCTION geo_coord_inside(this, poly) RESULT(inside)
862 TYPE(geo_coord), INTENT(IN) :: this
863 TYPE(geo_coordvect), INTENT(IN) :: poly
864 LOGICAL :: inside
865 
866 INTEGER :: i, j, starti
867 
868 inside = .false.
869 IF (all(poly%ll(1,:) == poly%ll(poly%vsize,:))) THEN ! Poligono chiuso
870  starti = 2
871  j = 1
872 ELSE ! Poligono non chiuso
873  starti = 1
874  j = poly%vsize
875 ENDIF
876 DO 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
888 ENDDO
889 
890 END FUNCTION geo_coord_inside
891 
892 
893 ELEMENTAL FUNCTION c_e_geo_coord(this) result (res)
894 TYPE(geo_coord),INTENT(in) :: this
895 LOGICAL :: res
896 
897 res = .not. this == geo_coord_miss
898 
899 end FUNCTION c_e_geo_coord
900 
901 
902 character(len=80) function to_char_geo_coord(this)
903 TYPE(geo_coord),INTENT(in) :: this
904 
905 to_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 
910 end function to_char_geo_coord
911 
912 
913 subroutine display_geo_coord(this)
914 TYPE(geo_coord),INTENT(in) :: this
915 
916 print*,trim(to_char(this))
918 end subroutine display_geo_coord
919 
920 
921 END MODULE geo_coord_class
Read a single geo_coord object or an array of geo_coord objects from a Fortran FORMATTED or UNFORMATT...
Constructors for the two classes.
Represent geo_coord object in a pretty string.
Derived type defining a one-dimensional array of georeferenced points with an associated topology (is...
Import one or more geo_coordvect objects from a plain text file or for a file in ESRI/Shapefile forma...
Utilities for managing files.
Export one or more geo_coordvect objects to a plain text file or to a file in ESRI/Shapefile format...
Classes for handling georeferenced sparse points in geographical corodinates.
Determine whether a point lies inside a polygon or a rectangle.
Gestione degli errori.
Definitions of constants and functions for working with missing values.
Costanti fisiche (DOUBLEPRECISION).
Definition: phys_const.f90:72
Write a single geo_coord object or an array of geo_coord objects to a Fortran FORMATTED or UNFORMATTE...
Derived type defining an isolated georeferenced point on Earth in polar geographical coordinates...
Utilities for CHARACTER variables.
Methods for returning the value of object members.
Detructors for the two classes.
Definition of constants to be used for declaring variables of a desired type.
Definition: kinds.F90:255

Generated with Doxygen.