50 INTEGER,
PARAMETER :: fp_geo=fp_d
51 real(kind=fp_geo),
PARAMETER :: fp_geo_miss=dmiss
57 INTEGER(kind=int_l) :: ilon, ilat
62 REAL(kind=fp_geo) :: lon, lat
71 REAL(kind=fp_geo),
POINTER :: ll(:,:) => null()
72 INTEGER :: vsize, vtype
76 TYPE(geo_coord
),
PARAMETER :: geo_coord_miss=geo_coord(imiss,imiss)
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
88 MODULE PROCEDURE geo_coord_init, geo_coordvect_init
95 MODULE PROCEDURE geo_coord_delete, geo_coordvect_delete
100 MODULE PROCEDURE geo_coord_getval, geo_coordvect_getval
104 INTERFACE operator (==)
105 MODULE PROCEDURE geo_coord_eq
109 INTERFACE operator (/=)
110 MODULE PROCEDURE geo_coord_ne
113 INTERFACE operator (>=)
114 MODULE PROCEDURE geo_coord_ge
117 INTERFACE operator (>)
118 MODULE PROCEDURE geo_coord_gt
121 INTERFACE operator (<=)
122 MODULE PROCEDURE geo_coord_le
125 INTERFACE operator (<)
126 MODULE PROCEDURE geo_coord_lt
132 MODULE PROCEDURE geo_coordvect_import, geo_coordvect_importvect
138 MODULE PROCEDURE geo_coordvect_export, geo_coordvect_exportvect
144 MODULE PROCEDURE geo_coord_read_unit, geo_coord_vect_read_unit
150 MODULE PROCEDURE geo_coord_write_unit, geo_coord_vect_write_unit
155 MODULE PROCEDURE geo_coord_inside, geo_coord_inside_rectang
160 MODULE PROCEDURE c_e_geo_coord
165 MODULE PROCEDURE to_char_geo_coord
170 MODULE PROCEDURE display_geo_coord
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
192 real(kind=fp_geo) :: llon,llat
194 CALL optio(ilon, this%ilon)
195 CALL optio(ilat, this%ilat)
197 if (.not.
c_e(this%ilon))
then
198 CALL optio(lon, llon)
199 if (
c_e(llon)) this%ilon=nint(llon*1.d5)
202 if (.not.
c_e(this%ilat))
then
203 CALL optio(lat, llat)
204 if (
c_e(llat)) this%ilat=nint(llat*1.d5)
207 END SUBROUTINE geo_coord_init
210 SUBROUTINE geo_coord_delete(this)
211 TYPE(geo_coord
),
INTENT(INOUT) :: this
216 END SUBROUTINE geo_coord_delete
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
229 IF (present(ilon)) ilon = getilon(this)
230 IF (present(ilat)) ilat = getilat(this)
232 IF (present(lon)) lon = getlon(this)
233 IF (present(lat)) lat = getlat(this)
235 END SUBROUTINE geo_coord_getval
242 elemental FUNCTION getilat(this)
243 TYPE(geo_coord
),
INTENT(IN) :: this
244 integer(kind=int_l) :: getilat
254 elemental FUNCTION getlat(this)
255 TYPE(geo_coord
),
INTENT(IN) :: this
256 real(kind=fp_geo) :: getlat
257 integer(kind=int_l) :: ilat
272 elemental FUNCTION getilon(this)
273 TYPE(geo_coord
),
INTENT(IN) :: this
274 integer(kind=int_l) :: getilon
285 elemental FUNCTION getlon(this)
286 TYPE(geo_coord
),
INTENT(IN) :: this
287 real(kind=fp_geo) :: getlon
288 integer(kind=int_l) :: ilon
300 elemental FUNCTION geo_coord_eq(this, that) RESULT(res)
301 TYPE(geo_coord
),
INTENT(IN) :: this, that
304 res = (this%ilon == that%ilon .AND. this%ilat == that%ilat)
306 END FUNCTION geo_coord_eq
309 elemental FUNCTION geo_coord_ne(this, that) RESULT(res)
310 TYPE(geo_coord
),
INTENT(IN) :: this, that
313 res = .not. this == that
315 END FUNCTION geo_coord_ne
319 elemental FUNCTION geo_coord_gt(this, that) RESULT(res)
320 TYPE(geo_coord
),
INTENT(IN) :: this, that
323 res = this%ilon > that%ilon
325 if ( this%ilon == that%ilon )
then
326 res = this%ilat > that%ilat
329 END FUNCTION geo_coord_gt
334 elemental FUNCTION geo_coord_ge(this, that) RESULT(res)
335 TYPE(geo_coord
),
INTENT(IN) :: this, that
338 res = .not. this < that
340 END FUNCTION geo_coord_ge
344 elemental FUNCTION geo_coord_lt(this, that) RESULT(res)
345 TYPE(geo_coord
),
INTENT(IN) :: this, that
348 res = this%ilon < that%ilon
350 if ( this%ilon == that%ilon )
then
351 res = this%ilat < that%ilat
355 END FUNCTION geo_coord_lt
360 elemental FUNCTION geo_coord_le(this, that) RESULT(res)
361 TYPE(geo_coord
),
INTENT(IN) :: this, that
364 res = .not. this > that
366 END FUNCTION geo_coord_le
371 elemental FUNCTION geo_coord_ure(this, that) RESULT(res)
372 TYPE(geo_coord
),
INTENT(IN) :: this, that
375 res = (this%ilon >= that%ilon .AND. this%ilat >= that%ilat)
377 END FUNCTION geo_coord_ure
381 elemental FUNCTION geo_coord_ur(this, that) RESULT(res)
382 TYPE(geo_coord
),
INTENT(IN) :: this, that
385 res = (this%ilon > that%ilon .AND. this%ilat > that%ilat)
387 END FUNCTION geo_coord_ur
392 elemental FUNCTION geo_coord_lle(this, that) RESULT(res)
393 TYPE(geo_coord
),
INTENT(IN) :: this, that
396 res = (this%ilon <= that%ilon .AND. this%ilat <= that%ilat)
398 END FUNCTION geo_coord_lle
402 elemental FUNCTION geo_coord_ll(this, that) RESULT(res)
403 TYPE(geo_coord
),
INTENT(IN) :: this, that
406 res = (this%ilon < that%ilon .AND. this%ilat < that%ilat)
408 END FUNCTION geo_coord_ll
416 SUBROUTINE geo_coord_read_unit(this, unit)
417 TYPE(geo_coord
),
INTENT(out) :: this
418 INTEGER,
INTENT(in) :: unit
420 CALL geo_coord_vect_read_unit((/this/), unit)
422 END SUBROUTINE geo_coord_read_unit
430 SUBROUTINE geo_coord_vect_read_unit(this, unit)
431 TYPE(geo_coord
) :: this(:)
432 INTEGER,
INTENT(in) :: unit
434 CHARACTER(len=40) :: form
437 INQUIRE(unit, form=form)
438 IF (form ==
'FORMATTED')
THEN
439 read(unit,*) (this(i)%ilon,this(i)%ilat, i=1,
SIZE(this))
443 READ(unit) (this(i)%ilon,this(i)%ilat, i=1,
SIZE(this))
446 END SUBROUTINE geo_coord_vect_read_unit
453 SUBROUTINE geo_coord_write_unit(this, unit)
454 TYPE(geo_coord
),
INTENT(in) :: this
455 INTEGER,
INTENT(in) :: unit
457 CALL geo_coord_vect_write_unit((/this/), unit)
459 END SUBROUTINE geo_coord_write_unit
466 SUBROUTINE geo_coord_vect_write_unit(this, unit)
467 TYPE(geo_coord
),
INTENT(in) :: this(:)
468 INTEGER,
INTENT(in) :: unit
470 CHARACTER(len=40) :: form
473 INQUIRE(unit, form=form)
474 IF (form ==
'FORMATTED')
THEN
475 WRITE(unit,*) (this(i)%ilon,this(i)%ilat, i=1,
SIZE(this))
479 WRITE(unit) (this(i)%ilon,this(i)%ilat, i=1,
SIZE(this))
482 END SUBROUTINE geo_coord_vect_write_unit
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
493 REAL(kind=fp_geo) :: x,y
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
500 END FUNCTION geo_coord_dist
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
517 res = (geo_coord_ure(this, coordmin) .AND. geo_coord_lle(this,coordmax))
519 END FUNCTION geo_coord_inside_rectang
533 RECURSIVE SUBROUTINE geo_coordvect_init(this, lon, lat)
535 REAL(kind=fp_geo),
INTENT(IN),
OPTIONAL :: lon(:)
536 REAL(kind=fp_geo),
INTENT(IN),
OPTIONAL :: lat(:)
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)
549 END SUBROUTINE geo_coordvect_init
554 SUBROUTINE geo_coordvect_delete(this)
557 IF (
ASSOCIATED(this%ll))
DEALLOCATE(this%ll)
561 END SUBROUTINE geo_coordvect_delete
572 SUBROUTINE geo_coordvect_getval(this, lon, lat)
574 REAL(kind=fp_geo),
OPTIONAL,
POINTER :: lon(:)
575 REAL(kind=fp_geo),
OPTIONAL,
POINTER :: lat(:)
577 IF (present(lon))
THEN
578 IF (
ASSOCIATED(this%ll))
THEN
579 ALLOCATE(lon(this%vsize))
580 lon(:) = this%ll(1:this%vsize,1)
583 IF (present(lat))
THEN
584 IF (
ASSOCIATED(this%ll))
THEN
585 ALLOCATE(lat(this%vsize))
586 lat(:) = this%ll(1:this%vsize,2)
590 END SUBROUTINE geo_coordvect_getval
593 SUBROUTINE geo_coordvect_import(this, unitsim, &
599 INTEGER,
OPTIONAL,
INTENT(IN) :: unitsim
601 TYPE(shpfileobject
),
OPTIONAL,
INTENT(INOUT) :: shphandle
603 INTEGER,
OPTIONAL,
INTENT(IN) :: nshp
605 REAL(kind=fp_geo),
ALLOCATABLE :: llon(:), llat(:)
606 REAL(kind=fp_geo) :: lv1,lv2,lv3,lv4
608 CHARACTER(len=40) :: lname
610 TYPE(shpobject
) :: shpobj
613 IF (present(unitsim))
THEN
615 READ(unitsim,*,end=10)lvsize,lv1,lv2,lv3,lv4,lname
616 ALLOCATE(llon(lvsize+1), llat(lvsize+1))
618 READ(unitsim,*)(llon(i),llat(i), i=1,lvsize)
620 IF (llon(1) /= llon(lvsize) .OR. llat(1) /= llat(lvsize))
THEN
622 llon(lvsize) = llon(1)
623 llat(lvsize) = llat(1)
626 CALL
init(this, lon=llon(1:lvsize), lat=llat(1:lvsize))
627 this%vtype = geo_coordvect_polygon
629 DEALLOCATE(llon, llat)
631 10 CALL raise_error(
'nella lettura del file '//trim(
to_char(unitsim)))
632 DEALLOCATE(llon, llat)
634 ELSE IF (present(shphandle) .AND. present(nshp))
THEN
636 shpobj = shpreadobject(shphandle, nshp)
637 IF (.NOT.shpisnull(shpobj))
THEN
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)
649 END SUBROUTINE geo_coordvect_import
652 SUBROUTINE geo_coordvect_export(this, unitsim, &
658 INTEGER,
OPTIONAL,
INTENT(IN) :: unitsim
660 TYPE(shpfileobject
),
OPTIONAL,
INTENT(INOUT) :: shphandle
662 INTEGER,
OPTIONAL,
INTENT(IN) :: nshp
666 TYPE(shpobject
) :: shpobj
669 IF (present(unitsim))
THEN
670 IF (this%vsize > 0)
THEN
672 WRITE(unitsim,*)
SIZE(this%ll,1),-1.,5000.,-0.1,1.1,
'Area'
674 WRITE(unitsim,*)(this%ll(i,1:2), i=1,this%vsize)
676 CALL raise_warning(
'oggetto geo_coordvect vuoto, non scrivo niente in '// &
680 ELSE IF (present(shphandle))
THEN
681 IF (present(nshp))
THEN
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
692 i=shpwriteobject(shphandle, lnshp, shpobj)
693 CALL shpdestroyobject(shpobj)
698 END SUBROUTINE geo_coordvect_export
718 SUBROUTINE geo_coordvect_importvect(this, shpfilesim, shpfile)
720 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: shpfilesim
721 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: shpfile
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
728 TYPE(shpfileobject
) :: shphandle
733 IF (present(shpfilesim))
THEN
735 OPEN(u, file=shpfilesim, status=
'old', err=30)
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)
747 CALL
import(this(i), unitsim=u)
752 IF (.NOT.
ASSOCIATED(this))
THEN
753 CALL raise_warning(
'file '//trim(shpfilesim)//
' vuoto o corrotto')
757 CALL raise_error(
'Impossibile aprire il file '//trim(shpfile))
760 ELSE IF (present(shpfile))
THEN
762 shphandle = shpopen(trim(shpfile),
'rb')
763 IF (shpfileisnull(shphandle))
THEN
764 CALL raise_error(
'Impossibile aprire lo shapefile '//trim(shpfile))
767 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
770 this(:)%vtype = shptype
772 CALL
import(this(i), shphandle=shphandle, nshp=i-1)
775 CALL shpclose(shphandle)
780 END SUBROUTINE geo_coordvect_importvect
785 SUBROUTINE geo_coordvect_exportvect(this, shpfilesim, shpfile, append)
787 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: shpfilesim
788 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: shpfile
794 LOGICAL,
INTENT(in),
OPTIONAL ::
append
796 REAL(kind=fp_d) :: minb(4), maxb(4)
797 INTEGER :: i, u, ns, shptype, dbfnf, dbfnr
800 TYPE(shpfileobject
) :: shphandle
808 IF (present(shpfilesim))
THEN
811 OPEN(u, file=shpfilesim, status=
'unknown', position=
'append', err=30)
813 OPEN(u, file=shpfilesim, status=
'unknown', err=30)
816 CALL
export(this(i), unitsim=u)
821 CALL raise_error(
'Impossibile aprire il file '//trim(shpfile))
823 ELSE IF (present(shpfile))
THEN
826 shphandle = shpopen(trim(shpfile),
'r+b')
827 IF (shpfileisnull(shphandle))
THEN
828 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
831 shphandle = shpcreate(trim(shpfile), geo_coordvect_polygon)
833 IF (shpfileisnull(shphandle))
THEN
834 CALL raise_error(
'Impossibile aprire lo shapefile '//trim(shpfile))
837 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
839 IF (i > ns .OR. lappend)
THEN
840 CALL
export(this(i), shphandle=shphandle)
842 CALL
export(this(i), shphandle=shphandle, nshp=i-1)
845 CALL shpclose(shphandle)
850 END SUBROUTINE geo_coordvect_exportvect
861 FUNCTION geo_coord_inside(this, poly) RESULT(inside)
862 TYPE(geo_coord
),
INTENT(IN) :: this
866 INTEGER :: i, j, starti
869 IF (all(poly%ll(1,:) == poly%ll(poly%vsize,:)))
THEN
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
890 END FUNCTION geo_coord_inside
893 ELEMENTAL FUNCTION c_e_geo_coord(this) result (res)
894 TYPE(geo_coord
),
INTENT(in) :: this
897 res = .not. this == geo_coord_miss
899 end FUNCTION c_e_geo_coord
902 character(len=80) function to_char_geo_coord(this)
903 TYPE(geo_coord
),
INTENT(in) :: this
905 to_char_geo_coord =
"GEO_COORD: Lon="// &
906 trim(
to_char(getlon(this),miss=
"Missing lon",form=
"(f11.5)"))//&
908 trim(
to_char(getlat(this),miss=
"Missing lat",form=
"(f11.5)"))
910 end function to_char_geo_coord
913 subroutine display_geo_coord(this)
914 TYPE(geo_coord
),
INTENT(in) :: this
918 end subroutine display_geo_coord
Definitions of constants and functions for working with missing values.
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...
Classes for handling georeferenced sparse points in geographical corodinates.
Export one or more geo_coordvect objects to a plain text file or to a file in ESRI/Shapefile format...
Determine whether a point lies inside a polygon or a rectangle.
Quick method to append an element to the array.
Utilities for CHARACTER variables.
Write a single geo_coord object or an array of geo_coord objects to a Fortran FORMATTED or UNFORMATTE...
Definition of constants to be used for declaring variables of a desired type.
Methods for returning the value of object members.
Utilities for managing files.
Detructors for the two classes.
Costanti fisiche (DOUBLEPRECISION).