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
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
183 CALL
optio(x, this%x)
184 CALL
optio(y, this%y)
186 END FUNCTION georef_coord_new
189 SUBROUTINE georef_coord_delete(this)
195 END SUBROUTINE georef_coord_delete
198 ELEMENTAL FUNCTION georef_coord_c_e(this) RESULT (res)
202 res = .NOT. this == georef_coord_miss
204 END FUNCTION georef_coord_c_e
213 ELEMENTAL SUBROUTINE georef_coord_getval(this, x, y)
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)
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)
256 DOUBLE PRECISION :: getlat
263 ELEMENTAL FUNCTION getlon(this)
265 DOUBLE PRECISION :: getlon
272 ELEMENTAL FUNCTION georef_coord_eq(this, that) RESULT(res)
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)
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)
303 res = .NOT.(this == that)
305 END FUNCTION georef_coord_ne
313 SUBROUTINE georef_coord_read_unit(this, unit)
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)
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)
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)
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)
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
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)
452 END SUBROUTINE georef_coord_array_delete
455 ELEMENTAL FUNCTION georef_coord_array_c_e(this) RESULT (res)
459 res =
ALLOCATED(this%coord)
461 END FUNCTION georef_coord_array_c_e
468 SUBROUTINE georef_coord_array_getval(this, x, y, topo, proj)
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
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.
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. ...