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)
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)
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)
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)
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)
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)
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, &
657 TYPE(geo_coordvect),
INTENT(INOUT) :: this
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)
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
796 REAL(kind=fp_d) :: minb(4), maxb(4)
797 INTEGER :: i, u, ns, shptype, dbfnf, dbfnr
800 TYPE(shpfileobject) :: shphandle
803 IF (
PRESENT(append))
THEN 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
863 TYPE(geo_coordvect),
INTENT(IN) :: poly
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
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.
Definitions of constants and functions for working with missing values.
Costanti fisiche (DOUBLEPRECISION).
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.