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)
274 integer(kind=int_l) :: getilon
285 elemental FUNCTION getlon(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)
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, &
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
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)
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)
897 res = .not. this == geo_coord_miss
899 end FUNCTION c_e_geo_coord
902 character(len=80) function to_char_geo_coord(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)
918 end 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...