49 DOUBLE PRECISION :: x=dmiss, y=dmiss
53 TYPE(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.
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
186 TYPE(georef_coord) :: this
188 CALL optio(x, this%x)
189 CALL optio(y, this%y)
191 END FUNCTION georef_coord_new
194 SUBROUTINE georef_coord_delete(this)
195 TYPE(georef_coord),
INTENT(inout) :: this
200 END SUBROUTINE georef_coord_delete
203 ELEMENTAL FUNCTION georef_coord_c_e(this)
RESULT (res)
204 TYPE(georef_coord),
INTENT(in) :: this
207 res = .NOT. this == georef_coord_miss
209 END FUNCTION georef_coord_c_e
218 ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
219 TYPE(georef_coord),
INTENT(in) :: this
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)
238 TYPE(georef_coord),
INTENT(in) :: this
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)
260 TYPE(georef_coord),
INTENT(in) :: this
261 DOUBLE PRECISION :: getlat
268 ELEMENTAL FUNCTION getlon(this)
270 DOUBLE PRECISION :: getlon
277 ELEMENTAL FUNCTION georef_coord_eq(this, that)
RESULT(res)
278 TYPE(georef_coord),
INTENT(IN) :: this, that
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)
287 TYPE(georef_coord),
INTENT(IN) :: this, that
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)
296 TYPE(georef_coord),
INTENT(IN) :: this, that
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)
305 TYPE(georef_coord),
INTENT(IN) :: this, that
308 res = .NOT.(this == that)
310 END FUNCTION georef_coord_ne
318 SUBROUTINE georef_coord_read_unit(this, unit)
319 TYPE(georef_coord),
INTENT(out) :: this
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)
333 TYPE(georef_coord) :: this(:)
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)
356 TYPE(georef_coord),
INTENT(in) :: this
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)
369 TYPE(georef_coord),
INTENT(in) :: this(:)
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)
392 TYPE(georef_coord),
INTENT (IN) :: that
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)
411 TYPE(georef_coord),
INTENT(IN) :: this
412 TYPE(georef_coord),
INTENT(IN) :: coordmin
413 TYPE(georef_coord),
INTENT(IN) :: coordmax
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
434 TYPE(georef_coord_array) :: this
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)
451 TYPE(georef_coord_array),
INTENT(inout) :: this
453 TYPE(georef_coord_array) :: lobj
457 END SUBROUTINE georef_coord_array_delete
460 ELEMENTAL FUNCTION georef_coord_array_c_e(this)
RESULT (res)
461 TYPE(georef_coord_array),
INTENT(in) :: this
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
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.