50INTEGER,
PARAMETER :: fp_geo=fp_d
51real(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
76TYPE(geo_coord),
PARAMETER :: geo_coord_miss=
geo_coord(imiss,imiss)
78INTEGER,
PARAMETER :: geo_coordvect_point = 1
79INTEGER,
PARAMETER :: geo_coordvect_arc = 3
80INTEGER,
PARAMETER :: geo_coordvect_polygon = 5
81INTEGER,
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
104INTERFACE OPERATOR (==)
105 MODULE PROCEDURE geo_coord_eq
109INTERFACE OPERATOR (/=)
110 MODULE PROCEDURE geo_coord_ne
113INTERFACE OPERATOR (>=)
114 MODULE PROCEDURE geo_coord_ge
117INTERFACE OPERATOR (>)
118 MODULE PROCEDURE geo_coord_gt
121INTERFACE OPERATOR (<=)
122 MODULE PROCEDURE geo_coord_le
125INTERFACE 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
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
192real(kind=fp_geo) :: llon,llat
194CALL optio(ilon, this%ilon)
195CALL optio(ilat, this%ilat)
197if (.not.
c_e(this%ilon))
then
198 CALL optio(lon, llon)
199 if (
c_e(llon)) this%ilon=nint(llon*1.d5)
202if (.not.
c_e(this%ilat))
then
203 CALL optio(lat, llat)
204 if (
c_e(llat)) this%ilat=nint(llat*1.d5)
207END SUBROUTINE geo_coord_init
210SUBROUTINE geo_coord_delete(this)
211TYPE(geo_coord),
INTENT(INOUT) :: this
216END SUBROUTINE geo_coord_delete
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
229IF (
PRESENT(ilon)) ilon = getilon(this)
230IF (
PRESENT(ilat)) ilat = getilat(this)
232IF (
PRESENT(lon)) lon = getlon(this)
233IF (
PRESENT(lat)) lat = getlat(this)
235END SUBROUTINE geo_coord_getval
242elemental FUNCTION getilat(this)
243TYPE(geo_coord),
INTENT(IN) :: this
244integer(kind=int_l) :: getilat
254elemental FUNCTION getlat(this)
255TYPE(geo_coord),
INTENT(IN) :: this
256real(kind=fp_geo) :: getlat
257integer(kind=int_l) :: ilat
272elemental FUNCTION getilon(this)
274integer(kind=int_l) :: getilon
285elemental FUNCTION getlon(this)
286TYPE(geo_coord),
INTENT(IN) :: this
287real(kind=fp_geo) :: getlon
288integer(kind=int_l) :: ilon
300elemental FUNCTION geo_coord_eq(this, that)
RESULT(res)
301TYPE(geo_coord),
INTENT(IN) :: this, that
304res = (this%ilon == that%ilon .AND. this%ilat == that%ilat)
306END FUNCTION geo_coord_eq
309elemental FUNCTION geo_coord_ne(this, that)
RESULT(res)
310TYPE(geo_coord),
INTENT(IN) :: this, that
313res = .not. this == that
315END FUNCTION geo_coord_ne
319elemental FUNCTION geo_coord_gt(this, that)
RESULT(res)
320TYPE(geo_coord),
INTENT(IN) :: this, that
323res = this%ilon > that%ilon
325if ( this%ilon == that%ilon )
then
326 res = this%ilat > that%ilat
329END FUNCTION geo_coord_gt
334elemental FUNCTION geo_coord_ge(this, that)
RESULT(res)
335TYPE(geo_coord),
INTENT(IN) :: this, that
338res = .not. this < that
340END FUNCTION geo_coord_ge
344elemental FUNCTION geo_coord_lt(this, that)
RESULT(res)
345TYPE(geo_coord),
INTENT(IN) :: this, that
348res = this%ilon < that%ilon
350if ( this%ilon == that%ilon )
then
351 res = this%ilat < that%ilat
355END FUNCTION geo_coord_lt
360elemental FUNCTION geo_coord_le(this, that)
RESULT(res)
361TYPE(geo_coord),
INTENT(IN) :: this, that
364res = .not. this > that
366END FUNCTION geo_coord_le
371elemental FUNCTION geo_coord_ure(this, that)
RESULT(res)
372TYPE(geo_coord),
INTENT(IN) :: this, that
375res = (this%ilon >= that%ilon .AND. this%ilat >= that%ilat)
377END FUNCTION geo_coord_ure
381elemental FUNCTION geo_coord_ur(this, that)
RESULT(res)
385res = (this%ilon > that%ilon .AND. this%ilat > that%ilat)
387END FUNCTION geo_coord_ur
392elemental FUNCTION geo_coord_lle(this, that)
RESULT(res)
396res = (this%ilon <= that%ilon .AND. this%ilat <= that%ilat)
398END FUNCTION geo_coord_lle
402elemental FUNCTION geo_coord_ll(this, that)
RESULT(res)
406res = (this%ilon < that%ilon .AND. this%ilat < that%ilat)
408END FUNCTION geo_coord_ll
416SUBROUTINE geo_coord_read_unit(this, unit)
418INTEGER,
INTENT(in) :: unit
420CALL geo_coord_vect_read_unit((/this/), unit)
422END SUBROUTINE geo_coord_read_unit
430SUBROUTINE geo_coord_vect_read_unit(this, unit)
432INTEGER,
INTENT(in) :: unit
434CHARACTER(len=40) :: form
437INQUIRE(unit, form=form)
438IF (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))
446END SUBROUTINE geo_coord_vect_read_unit
453SUBROUTINE geo_coord_write_unit(this, unit)
455INTEGER,
INTENT(in) :: unit
457CALL geo_coord_vect_write_unit((/this/), unit)
459END SUBROUTINE geo_coord_write_unit
466SUBROUTINE geo_coord_vect_write_unit(this, unit)
468INTEGER,
INTENT(in) :: unit
470CHARACTER(len=40) :: form
473INQUIRE(unit, form=form)
474IF (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))
482END SUBROUTINE geo_coord_vect_write_unit
487FUNCTION geo_coord_dist(this, that)
RESULT(dist)
491REAL(kind=fp_geo) :: dist
493REAL(kind=fp_geo) :: x,y
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
500END FUNCTION geo_coord_dist
511FUNCTION geo_coord_inside_rectang(this, coordmin, coordmax)
RESULT(res)
517res = (geo_coord_ure(this, coordmin) .AND. geo_coord_lle(this,coordmax))
519END FUNCTION geo_coord_inside_rectang
533RECURSIVE SUBROUTINE geo_coordvect_init(this, lon, lat)
535REAL(kind=fp_geo),
INTENT(IN),
OPTIONAL :: lon(:)
536REAL(kind=fp_geo),
INTENT(IN),
OPTIONAL :: lat(:)
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)
549END SUBROUTINE geo_coordvect_init
554SUBROUTINE geo_coordvect_delete(this)
557IF (
ASSOCIATED(this%ll))
DEALLOCATE(this%ll)
561END SUBROUTINE geo_coordvect_delete
572SUBROUTINE geo_coordvect_getval(this, lon, lat)
574REAL(kind=fp_geo),
OPTIONAL,
POINTER :: lon(:)
575REAL(kind=fp_geo),
OPTIONAL,
POINTER :: lat(:)
577IF (
PRESENT(lon))
THEN
578 IF (
ASSOCIATED(this%ll))
THEN
579 ALLOCATE(lon(this%vsize))
580 lon(:) = this%ll(1:this%vsize,1)
583IF (
PRESENT(lat))
THEN
584 IF (
ASSOCIATED(this%ll))
THEN
585 ALLOCATE(lat(this%vsize))
586 lat(:) = this%ll(1:this%vsize,2)
590END SUBROUTINE geo_coordvect_getval
593SUBROUTINE geo_coordvect_import(this, unitsim, &
599INTEGER,
OPTIONAL,
INTENT(IN) :: unitsim
601TYPE(shpfileobject),
OPTIONAL,
INTENT(INOUT) :: shphandle
603INTEGER,
OPTIONAL,
INTENT(IN) :: nshp
605REAL(kind=fp_geo),
ALLOCATABLE :: llon(:), llat(:)
606REAL(kind=fp_geo) :: lv1,lv2,lv3,lv4
608CHARACTER(len=40) :: lname
610TYPE(shpobject) :: shpobj
613IF (
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)
63110
CALL raise_error(
'nella lettura del file '//trim(
to_char(unitsim)))
632 DEALLOCATE(llon, llat)
634ELSE 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)
649END SUBROUTINE geo_coordvect_import
652SUBROUTINE geo_coordvect_export(this, unitsim, &
658INTEGER,
OPTIONAL,
INTENT(IN) :: unitsim
660TYPE(shpfileobject),
OPTIONAL,
INTENT(INOUT) :: shphandle
662INTEGER,
OPTIONAL,
INTENT(IN) :: nshp
666TYPE(shpobject) :: shpobj
669IF (
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 '// &
680ELSE 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)
698END SUBROUTINE geo_coordvect_export
718SUBROUTINE geo_coordvect_importvect(this, shpfilesim, shpfile)
720CHARACTER(len=*),
INTENT(in),
OPTIONAL :: shpfilesim
721CHARACTER(len=*),
INTENT(in),
OPTIONAL :: shpfile
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
728TYPE(shpfileobject) :: shphandle
733IF (
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))
760ELSE 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)
780END SUBROUTINE geo_coordvect_importvect
785SUBROUTINE geo_coordvect_exportvect(this, shpfilesim, shpfile, append)
787CHARACTER(len=*),
INTENT(in),
OPTIONAL :: shpfilesim
788CHARACTER(len=*),
INTENT(in),
OPTIONAL :: shpfile
794LOGICAL,
INTENT(in),
OPTIONAL :: append
796REAL(kind=fp_d) :: minb(4), maxb(4)
797INTEGER :: i, u, ns, shptype, dbfnf, dbfnr
800TYPE(shpfileobject) :: shphandle
803IF (
PRESENT(append))
THEN
808IF (
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))
823ELSE 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)
850END SUBROUTINE geo_coordvect_exportvect
861FUNCTION geo_coord_inside(this, poly)
RESULT(inside)
866INTEGER :: i, j, starti
869IF (all(poly%ll(1,:) == poly%ll(poly%vsize,:)))
THEN
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
890END FUNCTION geo_coord_inside
893ELEMENTAL FUNCTION c_e_geo_coord(this)
result (res)
897res = .not. this == geo_coord_miss
899end FUNCTION c_e_geo_coord
902character(len=80) function to_char_geo_coord(this)
905to_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)"))
910end function to_char_geo_coord
913subroutine display_geo_coord(this)
918end subroutine display_geo_coord
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).
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.
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...