615 READ(unitsim,*,
END=10)lvsize,lv1,lv2,lv3,lv4,lname
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)
672 WRITE(unitsim,*)
SIZE(this%ll,1),-1.,5000.,-0.1,1.1,
'Area'
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