49 DOUBLE PRECISION :: x=dmiss, y=dmiss
61 INTEGER,
ALLOCATABLE :: parts(:)
64 TYPE(geo_proj
) ::
proj
65 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
66 LOGICAL :: bbox_updated=.false.
69 INTEGER,
PARAMETER :: georef_coord_array_point = 1
70 INTEGER,
PARAMETER :: georef_coord_array_arc = 3
71 INTEGER,
PARAMETER :: georef_coord_array_polygon = 5
72 INTEGER,
PARAMETER :: georef_coord_array_multipoint = 8
79 MODULE PROCEDURE georef_coord_delete, georef_coord_array_delete
84 MODULE PROCEDURE georef_coord_c_e, georef_coord_array_c_e
89 MODULE PROCEDURE georef_coord_getval, georef_coord_proj_getval, georef_coord_array_getval
92 INTERFACE compute_bbox
93 MODULE PROCEDURE georef_coord_array_compute_bbox
97 INTERFACE operator (==)
98 MODULE PROCEDURE georef_coord_eq
102 INTERFACE operator (/=)
103 MODULE PROCEDURE georef_coord_ne
108 INTERFACE operator (>=)
109 MODULE PROCEDURE georef_coord_ge
114 INTERFACE operator (<=)
115 MODULE PROCEDURE georef_coord_le
122 MODULE PROCEDURE arrayof_georef_coord_array_import
128 MODULE PROCEDURE arrayof_georef_coord_array_export
135 MODULE PROCEDURE georef_coord_read_unit, georef_coord_vect_read_unit
141 MODULE PROCEDURE georef_coord_write_unit, georef_coord_vect_write_unit
146 MODULE PROCEDURE georef_coord_inside, georef_coord_inside_rectang
151 MODULE PROCEDURE georef_coord_dist
154 #define ARRAYOF_ORIGTYPE TYPE(georef_coord_array)
155 #define ARRAYOF_TYPE arrayof_georef_coord_array
157 #define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x)
158 #include "arrayof_pre.F90"
165 georef_coord_array_polygon, georef_coord_array_multipoint, &
166 delete,
c_e,
getval, compute_bbox, operator(==), operator(/=), operator(>=), operator(<=), &
171 georef_coord_new, georef_coord_array_new
175 #include "arrayof_post.F90"
183 FUNCTION georef_coord_new(x, y) RESULT(this)
184 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: x
185 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: y
188 CALL
optio(x, this%x)
189 CALL
optio(y, this%y)
191 END FUNCTION georef_coord_new
194 SUBROUTINE georef_coord_delete(this)
200 END SUBROUTINE georef_coord_delete
203 ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
207 res = .NOT. this == georef_coord_miss
209 END FUNCTION georef_coord_c_e
218 ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
220 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: x
221 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: y
223 IF (present(x)) x = this%x
224 IF (present(y)) y = this%y
226 END SUBROUTINE georef_coord_getval
237 ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
239 TYPE(geo_proj),
INTENT(in) :: proj
240 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: x
241 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: y
242 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: lon
243 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: lat
245 DOUBLE PRECISION :: llon, llat
247 IF (present(x)) x = this%x
248 IF (present(y)) y = this%y
249 IF (present(lon) .OR. present(lat))
THEN
251 IF (present(lon)) lon = llon
252 IF (present(lat)) lat = llat
255 END SUBROUTINE georef_coord_proj_getval
259 ELEMENTAL FUNCTION getlat(this)
261 DOUBLE PRECISION :: getlat
268 ELEMENTAL FUNCTION getlon(this)
270 DOUBLE PRECISION :: getlon
277 ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
281 res = (this%x == that%x .AND. this%y == that%y)
283 END FUNCTION georef_coord_eq
286 ELEMENTAL FUNCTION georef_coord_ge(this, that) RESULT(res)
290 res = (this%x >= that%x .AND. this%y >= that%y)
292 END FUNCTION georef_coord_ge
295 ELEMENTAL FUNCTION georef_coord_le(this, that) RESULT(res)
299 res = (this%x <= that%x .AND. this%y <= that%y)
301 END FUNCTION georef_coord_le
304 ELEMENTAL FUNCTION georef_coord_ne(this, that) RESULT(res)
308 res = .NOT.(this == that)
310 END FUNCTION georef_coord_ne
318 SUBROUTINE georef_coord_read_unit(this, unit)
320 INTEGER,
INTENT(in) :: unit
322 CALL georef_coord_vect_read_unit((/this/), unit)
324 END SUBROUTINE georef_coord_read_unit
332 SUBROUTINE georef_coord_vect_read_unit(this, unit)
334 INTEGER,
INTENT(in) :: unit
336 CHARACTER(len=40) :: form
339 INQUIRE(unit, form=form)
340 IF (form ==
'FORMATTED')
THEN
341 read(unit,*) (this(i)%x,this(i)%y, i=1,
SIZE(this))
345 READ(unit) (this(i)%x,this(i)%y, i=1,
SIZE(this))
348 END SUBROUTINE georef_coord_vect_read_unit
355 SUBROUTINE georef_coord_write_unit(this, unit)
357 INTEGER,
INTENT(in) :: unit
359 CALL georef_coord_vect_write_unit((/this/), unit)
361 END SUBROUTINE georef_coord_write_unit
368 SUBROUTINE georef_coord_vect_write_unit(this, unit)
370 INTEGER,
INTENT(in) :: unit
372 CHARACTER(len=40) :: form
375 INQUIRE(unit, form=form)
376 IF (form ==
'FORMATTED')
THEN
377 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,
SIZE(this))
381 WRITE(unit) (this(i)%x,this(i)%y, i=1,
SIZE(this))
384 END SUBROUTINE georef_coord_vect_write_unit
389 FUNCTION georef_coord_dist(this, that) RESULT(dist)
393 DOUBLE PRECISION ::
dist
395 DOUBLE PRECISION :: x,y
398 x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
400 dist = sqrt(x**2 + y**2)*degrad*rearth
402 END FUNCTION georef_coord_dist
410 FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax) RESULT(res)
416 res = (this >= coordmin .AND. this <= coordmax)
418 END FUNCTION georef_coord_inside_rectang
429 FUNCTION georef_coord_array_new(x, y, topo, proj) RESULT(this)
430 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: x(:)
431 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: y(:)
432 INTEGER,
INTENT(in),
OPTIONAL :: topo
433 type(geo_proj),
INTENT(in),
OPTIONAL ::
proj
438 IF (present(x) .AND. present(y))
THEN
439 lsize = min(
SIZE(x),
SIZE(y))
440 ALLOCATE(this%coord(lsize))
441 this%coord(1:lsize)%x = x(1:lsize)
442 this%coord(1:lsize)%y = y(1:lsize)
444 this%topo = optio_l(topo)
447 END FUNCTION georef_coord_array_new
450 SUBROUTINE georef_coord_array_delete(this)
457 END SUBROUTINE georef_coord_array_delete
460 ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
464 res =
ALLOCATED(this%coord)
466 END FUNCTION georef_coord_array_c_e
473 SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
475 DOUBLE PRECISION,
OPTIONAL,
ALLOCATABLE,
INTENT(out) :: x(:)
476 DOUBLE PRECISION,
OPTIONAL,
ALLOCATABLE,
INTENT(out) :: y(:)
478 INTEGER,
OPTIONAL,
INTENT(out) :: topo
479 type(geo_proj),
OPTIONAL,
INTENT(out) ::
proj
483 IF (
ALLOCATED(this%coord))
THEN
488 IF (
ALLOCATED(this%coord))
THEN
492 IF (present(topo)) topo = this%topo
495 END SUBROUTINE georef_coord_array_getval
503 SUBROUTINE georef_coord_array_compute_bbox(this)
506 IF (
ALLOCATED(this%coord))
THEN
507 this%bbox(1)%x = minval(this%coord(:)%x)
508 this%bbox(1)%y = minval(this%coord(:)%y)
509 this%bbox(2)%x = maxval(this%coord(:)%x)
510 this%bbox(2)%y = maxval(this%coord(:)%y)
511 this%bbox_updated = .true.
514 END SUBROUTINE georef_coord_array_compute_bbox
518 SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
520 TYPE(shpfileobject
),
INTENT(INOUT) :: shphandle
521 INTEGER,
INTENT(IN) :: nshp
523 TYPE(shpobject
) :: shpobj
526 shpobj = shpreadobject(shphandle, nshp)
527 IF (.NOT.shpisnull(shpobj))
THEN
529 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
530 topo=shpobj%nshptype)
531 IF (shpobj%nparts > 1 .AND.
ASSOCIATED(shpobj%panpartstart))
THEN
532 this%parts = shpobj%panpartstart(:)
533 ELSE IF (
ALLOCATED(this%parts))
THEN
534 DEALLOCATE(this%parts)
536 CALL shpdestroyobject(shpobj)
537 CALL compute_bbox(this)
541 END SUBROUTINE georef_coord_array_import
545 SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
547 TYPE(shpfileobject
),
INTENT(inout) :: shphandle
548 INTEGER,
INTENT(IN) :: nshp
551 TYPE(shpobject
) :: shpobj
553 IF (
ALLOCATED(this%coord))
THEN
554 IF (
ALLOCATED(this%parts))
THEN
555 shpobj = shpcreateobject(this%topo, -1,
SIZE(this%parts), this%parts, &
556 this%parts,
SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
558 shpobj = shpcreatesimpleobject(this%topo,
SIZE(this%coord), &
559 this%coord(:)%x, this%coord(:)%y)
565 IF (.NOT.shpisnull(shpobj))
THEN
566 i = shpwriteobject(shphandle, nshp, shpobj)
567 CALL shpdestroyobject(shpobj)
570 END SUBROUTINE georef_coord_array_export
583 SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
585 CHARACTER(len=*),
INTENT(in) :: shpfile
587 REAL(kind=fp_d) :: minb(4), maxb(4)
588 INTEGER :: i, ns, shptype, dbfnf, dbfnr
589 TYPE(shpfileobject
) :: shphandle
591 shphandle = shpopen(trim(shpfile),
'rb')
592 IF (shpfileisnull(shphandle))
THEN
599 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
601 CALL
insert(this, nelem=ns)
603 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
607 CALL shpclose(shphandle)
611 END SUBROUTINE arrayof_georef_coord_array_import
619 SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
621 CHARACTER(len=*),
INTENT(in) :: shpfile
624 TYPE(shpfileobject
) :: shphandle
626 IF (this%arraysize > 0)
THEN
627 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
629 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
631 IF (shpfileisnull(shphandle))
THEN
637 DO i = 1, this%arraysize
638 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
641 CALL shpclose(shphandle)
643 END SUBROUTINE arrayof_georef_coord_array_export
657 FUNCTION georef_coord_inside(this, poly) RESULT(inside)
665 IF (.NOT.
c_e(this))
RETURN
666 IF (.NOT.
ALLOCATED(poly%coord))
RETURN
668 IF (poly%bbox_updated)
THEN
669 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2)))
RETURN
672 IF (
ALLOCATED(poly%parts))
THEN
673 DO i = 1,
SIZE(poly%parts)-1
675 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
676 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
678 IF (
SIZE(poly%parts) > 0)
THEN
680 poly%coord(poly%parts(i)+1:)%x, &
681 poly%coord(poly%parts(i)+1:)%y)
685 IF (
SIZE(poly%coord) < 1)
RETURN
686 inside = pointinpoly(this%x, this%y, &
687 poly%coord(:)%x, poly%coord(:)%y)
692 FUNCTION pointinpoly(x, y, px, py)
693 DOUBLE PRECISION,
INTENT(in) :: x, y, px(:), py(:)
694 LOGICAL :: pointinpoly
696 INTEGER :: i, j, starti
698 pointinpoly = .false.
700 IF (px(1) == px(
SIZE(px)) .AND. py(1) == py(
SIZE(px)))
THEN
707 DO i = starti,
SIZE(px)
708 IF ((py(i) <= y .AND. y < py(j)) .OR. &
709 (py(j) <= y .AND. y < py(i)))
THEN
710 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i))
THEN
711 pointinpoly = .NOT. pointinpoly
717 END FUNCTION pointinpoly
719 END FUNCTION georef_coord_inside
Definitions of constants and functions for working with missing values.
Method for removing elements of the array at a desired position.
Compute forward coordinate transformation from geographical system to projected system.
Compute the distance in m between two points.
Method for packing the array object reducing at a minimum the memory occupation, without destroying i...
Read a single georef_coord object or an array of georef_coord objects from a Fortran FORMATTED or UNF...
Write a single georef_coord object or an array of georef_coord objects to a Fortran FORMATTED or UNFO...
Quick method to append an element to the array.
Derived type defining a dynamically extensible array of TYPE(georef_coord_array) elements.
Determine whether a point lies inside a polygon or a rectangle.
Generic subroutine for checking OPTIONAL parameters.
Import an array of georef_coord_array objects from a file in ESRI/Shapefile format.
Method for inserting elements of the array at a desired position.
Derived type defining a one-dimensional array of georeferenced points with an associated topology (is...
Compute backward coordinate transformation from projected system to geographical system.
Methods for returning the value of object members.
This module defines objects describing georeferenced sparse points possibly with topology and project...
Detructors for the two classes.
Derive type defining a single georeferenced point, either in geodetic or in projected coordinates...
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Export an array of georef_coord_array objects to a file in ESRI/Shapefile format. ...
Costanti fisiche (DOUBLEPRECISION).