49 DOUBLE PRECISION :: x=dmiss, y=dmiss
53TYPE(georef_coord),
PARAMETER :: georef_coord_miss=
georef_coord(dmiss,dmiss)
61 INTEGER,
ALLOCATABLE :: parts(:)
62 TYPE(georef_coord),
ALLOCATABLE :: coord(:)
64 TYPE(geo_proj) :: proj
65 TYPE(georef_coord) :: bbox(2)=(/georef_coord_miss, georef_coord_miss/)
66 LOGICAL :: bbox_updated=.false.
69INTEGER,
PARAMETER :: georef_coord_array_point = 1
70INTEGER,
PARAMETER :: georef_coord_array_arc = 3
71INTEGER,
PARAMETER :: georef_coord_array_polygon = 5
72INTEGER,
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
93 MODULE PROCEDURE georef_coord_array_compute_bbox
97INTERFACE OPERATOR (==)
98 MODULE PROCEDURE georef_coord_eq
102INTERFACE OPERATOR (/=)
103 MODULE PROCEDURE georef_coord_ne
108INTERFACE OPERATOR (>=)
109 MODULE PROCEDURE georef_coord_ge
114INTERFACE 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"
183FUNCTION georef_coord_new(x, y)
RESULT(this)
184DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: x
185DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: y
186TYPE(georef_coord) :: this
191END FUNCTION georef_coord_new
194SUBROUTINE georef_coord_delete(this)
195TYPE(georef_coord),
INTENT(inout) :: this
200END SUBROUTINE georef_coord_delete
203ELEMENTAL FUNCTION georef_coord_c_e(this)
RESULT (res)
204TYPE(georef_coord),
INTENT(in) :: this
207res = .NOT. this == georef_coord_miss
209END FUNCTION georef_coord_c_e
218ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
219TYPE(georef_coord),
INTENT(in) :: this
220DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: x
221DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: y
223IF (
PRESENT(x)) x = this%x
224IF (
PRESENT(y)) y = this%y
226END SUBROUTINE georef_coord_getval
237ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
238TYPE(georef_coord),
INTENT(in) :: this
239TYPE(geo_proj),
INTENT(in) :: proj
240DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: x
241DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: y
242DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: lon
243DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: lat
245DOUBLE PRECISION :: llon, llat
247IF (
PRESENT(x)) x = this%x
248IF (
PRESENT(y)) y = this%y
249IF (
PRESENT(lon) .OR.
present(lat))
THEN
251 IF (
PRESENT(lon)) lon = llon
252 IF (
PRESENT(lat)) lat = llat
255END SUBROUTINE georef_coord_proj_getval
259ELEMENTAL FUNCTION getlat(this)
261DOUBLE PRECISION :: getlat
268ELEMENTAL FUNCTION getlon(this)
269TYPE(georef_coord),
INTENT(in) :: this
270DOUBLE PRECISION :: getlon
277ELEMENTAL FUNCTION georef_coord_eq(this, that)
RESULT(res)
278TYPE(georef_coord),
INTENT(IN) :: this, that
281res = (this%x == that%x .AND. this%y == that%y)
283END FUNCTION georef_coord_eq
286ELEMENTAL FUNCTION georef_coord_ge(this, that)
RESULT(res)
287TYPE(georef_coord),
INTENT(IN) :: this, that
290res = (this%x >= that%x .AND. this%y >= that%y)
292END FUNCTION georef_coord_ge
295ELEMENTAL FUNCTION georef_coord_le(this, that)
RESULT(res)
296TYPE(georef_coord),
INTENT(IN) :: this, that
299res = (this%x <= that%x .AND. this%y <= that%y)
301END FUNCTION georef_coord_le
304ELEMENTAL FUNCTION georef_coord_ne(this, that)
RESULT(res)
305TYPE(georef_coord),
INTENT(IN) :: this, that
308res = .NOT.(this == that)
310END FUNCTION georef_coord_ne
318SUBROUTINE georef_coord_read_unit(this, unit)
319TYPE(georef_coord),
INTENT(out) :: this
320INTEGER,
INTENT(in) :: unit
322CALL georef_coord_vect_read_unit((/this/), unit)
324END SUBROUTINE georef_coord_read_unit
332SUBROUTINE georef_coord_vect_read_unit(this, unit)
333TYPE(georef_coord) :: this(:)
334INTEGER,
INTENT(in) :: unit
336CHARACTER(len=40) :: form
339INQUIRE(unit, form=form)
340IF (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))
348END SUBROUTINE georef_coord_vect_read_unit
355SUBROUTINE georef_coord_write_unit(this, unit)
356TYPE(georef_coord),
INTENT(in) :: this
357INTEGER,
INTENT(in) :: unit
359CALL georef_coord_vect_write_unit((/this/), unit)
361END SUBROUTINE georef_coord_write_unit
368SUBROUTINE georef_coord_vect_write_unit(this, unit)
369TYPE(georef_coord),
INTENT(in) :: this(:)
370INTEGER,
INTENT(in) :: unit
372CHARACTER(len=40) :: form
375INQUIRE(unit, form=form)
376IF (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))
384END SUBROUTINE georef_coord_vect_write_unit
389FUNCTION georef_coord_dist(this, that)
RESULT(dist)
391TYPE(georef_coord),
INTENT (IN) :: this
392TYPE(georef_coord),
INTENT (IN) :: that
393DOUBLE PRECISION :: dist
395DOUBLE PRECISION :: x,y
398x = (this%x-that%x)*cos(((this%y+this%y)/2.)*degrad)
400dist = sqrt(x**2 + y**2)*degrad*rearth
402END FUNCTION georef_coord_dist
410FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax)
RESULT(res)
411TYPE(georef_coord),
INTENT(IN) :: this
412TYPE(georef_coord),
INTENT(IN) :: coordmin
413TYPE(georef_coord),
INTENT(IN) :: coordmax
416res = (this >= coordmin .AND. this <= coordmax)
418END FUNCTION georef_coord_inside_rectang
429FUNCTION georef_coord_array_new(x, y, topo, proj)
RESULT(this)
430DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: x(:)
431DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: y(:)
432INTEGER,
INTENT(in),
OPTIONAL :: topo
433TYPE(geo_proj),
INTENT(in),
OPTIONAL :: proj
434TYPE(georef_coord_array) :: this
438IF (
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)
444this%topo = optio_l(topo)
447END FUNCTION georef_coord_array_new
450SUBROUTINE georef_coord_array_delete(this)
451TYPE(georef_coord_array),
INTENT(inout) :: this
453TYPE(georef_coord_array) :: lobj
457END SUBROUTINE georef_coord_array_delete
460ELEMENTAL FUNCTION georef_coord_array_c_e(this)
RESULT (res)
461TYPE(georef_coord_array),
INTENT(in) :: this
464res =
ALLOCATED(this%coord)
466END FUNCTION georef_coord_array_c_e
473SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
475DOUBLE PRECISION,
OPTIONAL,
ALLOCATABLE,
INTENT(out) :: x(:)
476DOUBLE PRECISION,
OPTIONAL,
ALLOCATABLE,
INTENT(out) :: y(:)
478INTEGER,
OPTIONAL,
INTENT(out) :: topo
479TYPE(geo_proj),
OPTIONAL,
INTENT(out) :: proj
483 IF (
ALLOCATED(this%coord))
THEN
488 IF (
ALLOCATED(this%coord))
THEN
492IF (
PRESENT(topo)) topo = this%topo
495END SUBROUTINE georef_coord_array_getval
503SUBROUTINE georef_coord_array_compute_bbox(this)
506IF (
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.
514END SUBROUTINE georef_coord_array_compute_bbox
518SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
520TYPE(shpfileobject),
INTENT(INOUT) :: shphandle
521INTEGER,
INTENT(IN) :: nshp
523TYPE(shpobject) :: shpobj
526shpobj = shpreadobject(shphandle, nshp)
527IF (.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)
541END SUBROUTINE georef_coord_array_import
545SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
547TYPE(shpfileobject),
INTENT(inout) :: shphandle
548INTEGER,
INTENT(IN) :: nshp
551TYPE(shpobject) :: shpobj
553IF (
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)
565IF (.NOT.shpisnull(shpobj))
THEN
566 i = shpwriteobject(shphandle, nshp, shpobj)
567 CALL shpdestroyobject(shpobj)
570END SUBROUTINE georef_coord_array_export
583SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
585CHARACTER(len=*),
INTENT(in) :: shpfile
587REAL(kind=fp_d) :: minb(4), maxb(4)
588INTEGER :: i, ns, shptype, dbfnf, dbfnr
589TYPE(shpfileobject) :: shphandle
591shphandle = shpopen(trim(shpfile),
'rb')
592IF (shpfileisnull(shphandle))
THEN
599CALL 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)
607CALL shpclose(shphandle)
611END SUBROUTINE arrayof_georef_coord_array_import
619SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
621CHARACTER(len=*),
INTENT(in) :: shpfile
624TYPE(shpfileobject) :: shphandle
626IF (this%arraysize > 0)
THEN
627 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
629 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
631IF (shpfileisnull(shphandle))
THEN
637DO i = 1, this%arraysize
638 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
641CALL shpclose(shphandle)
643END SUBROUTINE arrayof_georef_coord_array_export
657FUNCTION georef_coord_inside(this, poly)
RESULT(inside)
665IF (.NOT.
c_e(this))
RETURN
666IF (.NOT.
ALLOCATED(poly%coord))
RETURN
668IF (poly%bbox_updated)
THEN
669 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2)))
RETURN
672IF (
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)
692FUNCTION pointinpoly(x, y, px, py)
693DOUBLE PRECISION,
INTENT(in) :: x, y, px(:), py(:)
694LOGICAL :: pointinpoly
696INTEGER :: i, j, starti
700IF (px(1) == px(
SIZE(px)) .AND. py(1) == py(
SIZE(px)))
THEN
707DO 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
717END FUNCTION pointinpoly
719END FUNCTION georef_coord_inside
Compute forward coordinate transformation from geographical system to projected system.
Compute backward coordinate transformation from projected system to geographical system.
Quick method to append an element to the array.
Detructors for the two classes.
Compute the distance in m between two points.
Export an array of georef_coord_array objects to a file in ESRI/Shapefile format.
Methods for returning the value of object members.
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.
Determine whether a point lies inside a polygon or a rectangle.
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...
Method for removing elements of the array at a desired position.
Write a single georef_coord object or an array of georef_coord objects to a Fortran FORMATTED or UNFO...
Generic subroutine for checking OPTIONAL parameters.
Costanti fisiche (DOUBLEPRECISION).
This module defines objects describing georeferenced sparse points possibly with topology and project...
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Derived type defining a dynamically extensible array of TYPE(georef_coord_array) elements.
Derived type defining a one-dimensional array of georeferenced points with an associated topology (is...
Derive type defining a single georeferenced point, either in geodetic or in projected coordinates.