libsim  Versione 7.1.7
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
267 
272 elemental FUNCTION getilon(this)
273 TYPE(geo_coord),INTENT(IN) :: this
274 integer(kind=int_l) :: getilon
275 
276 getilon = this%ilon
277 
278 END FUNCTION getilon
279 
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
307 
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
341 
344 elemental FUNCTION geo_coord_lt(this, that) RESULT(res)
345 TYPE(geo_coord),INTENT(IN) :: this, that
346 LOGICAL :: res
347 
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 
357 
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
367 
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
384 
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
409 
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)
421 
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 
484 
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))
518 
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
591 
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)
732 
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))
917 
918 end subroutine display_geo_coord
919 
920 
921 END 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.