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
149 #define ARRAYOF_ORIGTYPE TYPE(georef_coord_array) 150 #define ARRAYOF_TYPE arrayof_georef_coord_array 152 #define ARRAYOF_ORIGDESTRUCTOR(x) CALL delete(x) 153 #include "arrayof_pre.F90" 160 georef_coord_array_polygon, georef_coord_array_multipoint, &
161 delete,
c_e,
getval, compute_bbox,
OPERATOR(==),
OPERATOR(/=),
OPERATOR(>=),
OPERATOR(<=), &
166 georef_coord_new, georef_coord_array_new
170 #include "arrayof_post.F90" 178 FUNCTION georef_coord_new(x, y)
RESULT(this)
179 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: x
180 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: y
181 TYPE(georef_coord) :: this
183 CALL optio(x, this%x)
184 CALL optio(y, this%y)
186 END FUNCTION georef_coord_new
189 SUBROUTINE georef_coord_delete(this)
190 TYPE(georef_coord),
INTENT(inout) :: this
195 END SUBROUTINE georef_coord_delete
198 ELEMENTAL FUNCTION georef_coord_c_e(this)
RESULT (res)
199 TYPE(georef_coord),
INTENT(in) :: this
202 res = .NOT. this == georef_coord_miss
204 END FUNCTION georef_coord_c_e
213 ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
214 TYPE(georef_coord),
INTENT(in) :: this
215 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: x
216 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: y
218 IF (
PRESENT(x)) x = this%x
219 IF (
PRESENT(y)) y = this%y
221 END SUBROUTINE georef_coord_getval
232 ELEMENTAL SUBROUTINE georef_coord_proj_getval(this, proj, x, y, lon, lat)
233 TYPE(georef_coord),
INTENT(in) :: this
234 TYPE(geo_proj),
INTENT(in) :: proj
235 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: x
236 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: y
237 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: lon
238 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: lat
240 DOUBLE PRECISION :: llon, llat
242 IF (
PRESENT(x)) x = this%x
243 IF (
PRESENT(y)) y = this%y
244 IF (
PRESENT(lon) .OR.
present(lat))
THEN 246 IF (
PRESENT(lon)) lon = llon
247 IF (
PRESENT(lat)) lat = llat
250 END SUBROUTINE georef_coord_proj_getval
254 ELEMENTAL FUNCTION getlat(this)
255 TYPE(georef_coord),
INTENT(in) :: this
256 DOUBLE PRECISION :: getlat
263 ELEMENTAL FUNCTION getlon(this)
264 TYPE(georef_coord),
INTENT(in) :: this
265 DOUBLE PRECISION :: getlon
272 ELEMENTAL FUNCTION georef_coord_eq(this, that)
RESULT(res)
273 TYPE(georef_coord),
INTENT(IN) :: this, that
276 res = (this%x == that%x .AND. this%y == that%y)
278 END FUNCTION georef_coord_eq
281 ELEMENTAL FUNCTION georef_coord_ge(this, that)
RESULT(res)
285 res = (this%x >= that%x .AND. this%y >= that%y)
287 END FUNCTION georef_coord_ge
290 ELEMENTAL FUNCTION georef_coord_le(this, that)
RESULT(res)
291 TYPE(georef_coord),
INTENT(IN) :: this, that
294 res = (this%x <= that%x .AND. this%y <= that%y)
296 END FUNCTION georef_coord_le
299 ELEMENTAL FUNCTION georef_coord_ne(this, that)
RESULT(res)
300 TYPE(georef_coord),
INTENT(IN) :: this, that
303 res = .NOT.(this == that)
305 END FUNCTION georef_coord_ne
313 SUBROUTINE georef_coord_read_unit(this, unit)
314 TYPE(georef_coord),
INTENT(out) :: this
315 INTEGER,
INTENT(in) :: unit
317 CALL georef_coord_vect_read_unit((/this/), unit)
319 END SUBROUTINE georef_coord_read_unit
327 SUBROUTINE georef_coord_vect_read_unit(this, unit)
328 TYPE(georef_coord) :: this(:)
329 INTEGER,
INTENT(in) :: unit
331 CHARACTER(len=40) :: form
334 INQUIRE(unit, form=form)
335 IF (form ==
'FORMATTED')
THEN 336 read(unit,*) (this(i)%x,this(i)%y, i=1,
SIZE(this))
340 READ(unit) (this(i)%x,this(i)%y, i=1,
SIZE(this))
343 END SUBROUTINE georef_coord_vect_read_unit
350 SUBROUTINE georef_coord_write_unit(this, unit)
351 TYPE(georef_coord),
INTENT(in) :: this
352 INTEGER,
INTENT(in) :: unit
354 CALL georef_coord_vect_write_unit((/this/), unit)
356 END SUBROUTINE georef_coord_write_unit
363 SUBROUTINE georef_coord_vect_write_unit(this, unit)
364 TYPE(georef_coord),
INTENT(in) :: this(:)
365 INTEGER,
INTENT(in) :: unit
367 CHARACTER(len=40) :: form
370 INQUIRE(unit, form=form)
371 IF (form ==
'FORMATTED')
THEN 372 WRITE(unit,*) (this(i)%x,this(i)%y, i=1,
SIZE(this))
376 WRITE(unit) (this(i)%x,this(i)%y, i=1,
SIZE(this))
379 END SUBROUTINE georef_coord_vect_write_unit
405 FUNCTION georef_coord_inside_rectang(this, coordmin, coordmax)
RESULT(res)
406 TYPE(georef_coord),
INTENT(IN) :: this
407 TYPE(georef_coord),
INTENT(IN) :: coordmin
408 TYPE(georef_coord),
INTENT(IN) :: coordmax
411 res = (this >= coordmin .AND. this <= coordmax)
413 END FUNCTION georef_coord_inside_rectang
424 FUNCTION georef_coord_array_new(x, y, topo, proj)
RESULT(this)
425 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: x(:)
426 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: y(:)
427 INTEGER,
INTENT(in),
OPTIONAL :: topo
428 TYPE(geo_proj),
INTENT(in),
OPTIONAL :: proj
429 TYPE(georef_coord_array) :: this
433 IF (
PRESENT(x) .AND.
PRESENT(y))
THEN 434 lsize = min(
SIZE(x),
SIZE(y))
435 ALLOCATE(this%coord(lsize))
436 this%coord(1:lsize)%x = x(1:lsize)
437 this%coord(1:lsize)%y = y(1:lsize)
439 this%topo = optio_l(topo)
442 END FUNCTION georef_coord_array_new
445 SUBROUTINE georef_coord_array_delete(this)
446 TYPE(georef_coord_array),
INTENT(inout) :: this
448 TYPE(georef_coord_array) :: lobj
452 END SUBROUTINE georef_coord_array_delete
455 ELEMENTAL FUNCTION georef_coord_array_c_e(this)
RESULT (res)
456 TYPE(georef_coord_array),
INTENT(in) :: this
459 res =
ALLOCATED(this%coord)
461 END FUNCTION georef_coord_array_c_e
468 SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
469 TYPE(georef_coord_array),
INTENT(in) :: this
470 DOUBLE PRECISION,
OPTIONAL,
ALLOCATABLE,
INTENT(out) :: x(:)
471 DOUBLE PRECISION,
OPTIONAL,
ALLOCATABLE,
INTENT(out) :: y(:)
473 INTEGER,
OPTIONAL,
INTENT(out) :: topo
474 TYPE(geo_proj),
OPTIONAL,
INTENT(out) :: proj
478 IF (
ALLOCATED(this%coord))
THEN 483 IF (
ALLOCATED(this%coord))
THEN 487 IF (
PRESENT(topo)) topo = this%topo
490 END SUBROUTINE georef_coord_array_getval
498 SUBROUTINE georef_coord_array_compute_bbox(this)
501 IF (
ALLOCATED(this%coord))
THEN 502 this%bbox(1)%x = minval(this%coord(:)%x)
503 this%bbox(1)%y = minval(this%coord(:)%y)
504 this%bbox(2)%x = maxval(this%coord(:)%x)
505 this%bbox(2)%y = maxval(this%coord(:)%y)
506 this%bbox_updated = .true.
509 END SUBROUTINE georef_coord_array_compute_bbox
513 SUBROUTINE georef_coord_array_import(this, shphandle, nshp)
515 TYPE(shpfileobject),
INTENT(INOUT) :: shphandle
516 INTEGER,
INTENT(IN) :: nshp
518 TYPE(shpobject) :: shpobj
521 shpobj = shpreadobject(shphandle, nshp)
522 IF (.NOT.shpisnull(shpobj))
THEN 524 this = georef_coord_array_new(x=dble(shpobj%padfx), y=dble(shpobj%padfy), &
525 topo=shpobj%nshptype)
526 IF (shpobj%nparts > 1 .AND.
ASSOCIATED(shpobj%panpartstart))
THEN 527 this%parts = shpobj%panpartstart(:)
528 ELSE IF (
ALLOCATED(this%parts))
THEN 529 DEALLOCATE(this%parts)
531 CALL shpdestroyobject(shpobj)
532 CALL compute_bbox(this)
536 END SUBROUTINE georef_coord_array_import
540 SUBROUTINE georef_coord_array_export(this, shphandle, nshp)
542 TYPE(shpfileobject),
INTENT(inout) :: shphandle
543 INTEGER,
INTENT(IN) :: nshp
546 TYPE(shpobject) :: shpobj
548 IF (
ALLOCATED(this%coord))
THEN 549 IF (
ALLOCATED(this%parts))
THEN 550 shpobj = shpcreateobject(this%topo, -1,
SIZE(this%parts), this%parts, &
551 this%parts,
SIZE(this%coord), this%coord(:)%x, this%coord(:)%y)
553 shpobj = shpcreatesimpleobject(this%topo,
SIZE(this%coord), &
554 this%coord(:)%x, this%coord(:)%y)
560 IF (.NOT.shpisnull(shpobj))
THEN 561 i = shpwriteobject(shphandle, nshp, shpobj)
562 CALL shpdestroyobject(shpobj)
565 END SUBROUTINE georef_coord_array_export
578 SUBROUTINE arrayof_georef_coord_array_import(this, shpfile)
580 CHARACTER(len=*),
INTENT(in) :: shpfile
582 REAL(kind=fp_d) :: minb(4), maxb(4)
583 INTEGER :: i, ns, shptype, dbfnf, dbfnr
584 TYPE(shpfileobject) :: shphandle
586 shphandle = shpopen(trim(shpfile),
'rb')
587 IF (shpfileisnull(shphandle))
THEN 594 CALL shpgetinfo(shphandle, ns, shptype, minb, maxb, dbfnf, dbfnr)
596 CALL insert(this, nelem=ns)
598 CALL georef_coord_array_import(this%array(i), shphandle=shphandle, nshp=i-1)
602 CALL shpclose(shphandle)
606 END SUBROUTINE arrayof_georef_coord_array_import
614 SUBROUTINE arrayof_georef_coord_array_export(this, shpfile)
616 CHARACTER(len=*),
INTENT(in) :: shpfile
619 TYPE(shpfileobject) :: shphandle
621 IF (this%arraysize > 0)
THEN 622 shphandle = shpcreate(trim(shpfile), this%array(1)%topo)
624 shphandle = shpcreate(trim(shpfile), georef_coord_array_polygon)
626 IF (shpfileisnull(shphandle))
THEN 632 DO i = 1, this%arraysize
633 CALL georef_coord_array_export(this%array(i), shphandle=shphandle, nshp=i-1)
636 CALL shpclose(shphandle)
638 END SUBROUTINE arrayof_georef_coord_array_export
652 FUNCTION georef_coord_inside(this, poly)
RESULT(inside)
660 IF (.NOT.
c_e(this))
RETURN 661 IF (.NOT.
ALLOCATED(poly%coord))
RETURN 663 IF (poly%bbox_updated)
THEN 664 IF (.NOT.georef_coord_inside_rectang(this, poly%bbox(1), poly%bbox(2)))
RETURN 667 IF (
ALLOCATED(poly%parts))
THEN 668 DO i = 1,
SIZE(poly%parts)-1
670 poly%coord(poly%parts(i)+1:poly%parts(i+1))%x, &
671 poly%coord(poly%parts(i)+1:poly%parts(i+1))%y)
673 IF (
SIZE(poly%parts) > 0)
THEN 675 poly%coord(poly%parts(i)+1:)%x, &
676 poly%coord(poly%parts(i)+1:)%y)
680 IF (
SIZE(poly%coord) < 1)
RETURN 681 inside = pointinpoly(this%x, this%y, &
682 poly%coord(:)%x, poly%coord(:)%y)
687 FUNCTION pointinpoly(x, y, px, py)
688 DOUBLE PRECISION,
INTENT(in) :: x, y, px(:), py(:)
689 LOGICAL :: pointinpoly
691 INTEGER :: i, j, starti
693 pointinpoly = .false.
695 IF (px(1) == px(
SIZE(px)) .AND. py(1) == py(
SIZE(px)))
THEN 702 DO i = starti,
SIZE(px)
703 IF ((py(i) <= y .AND. y < py(j)) .OR. &
704 (py(j) <= y .AND. y < py(i)))
THEN 705 IF (x < (px(j) - px(i)) * (y - py(i)) / (py(j) - py(i)) + px(i))
THEN 706 pointinpoly = .NOT. pointinpoly
712 END FUNCTION pointinpoly
714 END FUNCTION georef_coord_inside
Method for removing elements of the array at a desired position.
Compute forward coordinate transformation from geographical system to projected system.
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...
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
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.
This module defines objects describing georeferenced sparse points possibly with topology and project...
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.
Definitions of constants and functions for working with missing values.
Detructors for the two classes.
Derive type defining a single georeferenced point, either in geodetic or in projected coordinates...
Export an array of georef_coord_array objects to a file in ESRI/Shapefile format. ...