FortranGIS Version 3.0
proj.F90
1! Copyright 2011 Davide Cesari <dcesari69 at gmail dot com>
2!
3! This file is part of FortranGIS.
4!
5! FortranGIS is free software: you can redistribute it and/or modify
6! it under the terms of the GNU Lesser General Public License as
7! published by the Free Software Foundation, either version 3 of the
8! License, or (at your option) any later version.
9!
10! FortranGIS is distributed in the hope that it will be useful, but
11! WITHOUT ANY WARRANTY; without even the implied warranty of
12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13! Lesser General Public License for more details.
14!
15! You should have received a copy of the GNU Lesser General Public
16! License along with FortranGIS. If not, see
17! <http://www.gnu.org/licenses/>.
18
19!> Fortran 2003 interface to the proj.4 https://github.com/OSGeo/proj.4 library.
20!! The following functions are directly interfaced to their
21!! corresponding C version, so they are undocumented here, please
22!! refer to the original proj C API documentation, e.g. at the address
23!! https://github.com/OSGeo/proj.4/wiki/ProjAPI , for their use:
24!! - pj_init_plus() -> FUNCTION pj_init_plus()
25!! - pj_transform()() -> FUNCTION pj_transform()()
26!! - pj_datum_transform() -> FUNCTION pj_datum_transform()
27!! - pj_fwd() -> FUNCTION pj_fwd()
28!! - pj_inv() -> FUNCTION pj_inv()
29!! - pj_geocentric_to_geodetic() -> FUNCTION pj_geocentric_to_geodetic()
30!! - pj_geodetic_to_geocentric() -> FUNCTION pj_geodetic_to_geocentric()
31!! - pj_compare_datums() -> FUNCTION pj_compare_datums()
32!! - pj_is_latlong() -> FUNCTION pj_is_latlong()
33!! - pj_is_geocent() -> FUNCTION pj_is_geocent()
34!! - pj_get_def() -> FUNCTION pj_get_def()
35!! - pj_latlong_from_proj() -> FUNCTION pj_latlong_from_proj()
36!! - pj_free() -> SUBROUTINE pj_free()
37!!
38!! Notice that, if relevant, the result of functions returning an
39!! integer has to be interpreted as 0=false, nonzero=true or 0=ok,
40!! nonzero=error.
41!!
42!! Some of these functions have also a more Fortran-friendly interface
43!! explicitely documented here, with an \a _f appended to the name.
44!!
45!! For an example of application of the \a proj module, please refer
46!! to the following test program, which performs a forward and
47!! backward transformation:
48!! \include proj_test.F90
49!!
50!! \ingroup libfortrangis
51MODULE proj
52use,INTRINSIC :: iso_c_binding
53IMPLICIT NONE
54
55!> Object describing a cartographic projection.
56!! Its components are private so they should not be manipulated
57!! directly.
58TYPE,BIND(C) :: pj_object
59 PRIVATE
60 TYPE(c_ptr) :: ptr = c_null_ptr
61END TYPE pj_object
63!> Object representing a null cartographic projection.
64!! It should be used for comparison of function results
65!! with the \a == operator for error checking.
66TYPE(pj_object),PARAMETER :: pj_object_null=pj_object(c_null_ptr)
67
68!> Object describing a coordinate pair.
69!! It can indicate either projected coordinate (u=x, v=y) or
70!! spherical coordinates (u=lon, v=lat).
71TYPE,BIND(C) :: pjuv_object
72 REAL(kind=c_double) :: u=huge(1.0_c_double)
73 REAL(kind=c_double) :: v=huge(1.0_c_double)
74END TYPE pjuv_object
75
76REAL(kind=c_double),PARAMETER :: pj_deg_to_rad=.0174532925199432958d0 !< equivalent to the C api symbol DEG_TO_RAD
77REAL(kind=c_double),PARAMETER :: pj_rad_to_deg=57.29577951308232d0 !< equivalent to the C api symbol RAD_TO_DEG
78
79!> Test whether an opaque object is valid.
80!! Please use the interface name pj_associated, but for the documentation
81!! see the specific function pj_associated_object.
82INTERFACE pj_associated
83 MODULE PROCEDURE pj_associated_object
84END INTERFACE pj_associated
85
86!> Initialize a projection from a string.
87!! It returns an object of pj_object type.
88INTERFACE
89 FUNCTION pj_init_plus(name) bind(C,name='pj_init_plus')
90 IMPORT
91 CHARACTER(kind=c_char) :: name(*) !< Projection string, must be terminated by //CHAR(0)
92 TYPE(pj_object) :: pj_init_plus
93 END FUNCTION pj_init_plus
94END INTERFACE
95
96INTERFACE
97 FUNCTION pj_transform(src, dst, point_count, point_offset, x, y, z) &
98 bind(c,name='pj_transform')
99 IMPORT
100 TYPE(pj_object),VALUE :: src
101 TYPE(pj_object),VALUE :: dst
102 INTEGER(kind=c_long),VALUE :: point_count
103 INTEGER(kind=c_int),VALUE :: point_offset
104 TYPE(c_ptr), VALUE :: x !< Must be C-pointer to an array real(kind=c_double) :: x(:)
105 TYPE(c_ptr), VALUE :: y !< Must be C-pointer to an array real(kind=c_double) :: y(:)
106 TYPE(c_ptr), VALUE :: z !< Must be C-pointer to an array real(kind=c_double) :: z(:)
107 INTEGER(kind=c_int) :: pj_transform
108 END FUNCTION pj_transform
109END INTERFACE
110
111INTERFACE
112 FUNCTION pj_datum_transform(src, dst, point_count, point_offset, x, y, z) &
113 bind(c,name='pj_datum_transform')
114 IMPORT
115 TYPE(pj_object),VALUE :: src
116 TYPE(pj_object),VALUE :: dst
117 INTEGER(kind=c_long),VALUE :: point_count
118 INTEGER(kind=c_int),VALUE :: point_offset
119 REAL(kind=c_double) :: x(*)
120 REAL(kind=c_double) :: y(*)
121 REAL(kind=c_double) :: z(*)
122 INTEGER(kind=c_int) :: pj_datum_transform
123 END FUNCTION pj_datum_transform
124END INTERFACE
125
126INTERFACE
127 FUNCTION pj_fwd(val, proj) bind(C,name='pj_fwd')
128 IMPORT
129 TYPE(pjuv_object),VALUE :: val
130 TYPE(pj_object),VALUE :: proj
131 TYPE(pjuv_object) :: pj_fwd
132 END FUNCTION pj_fwd
133END INTERFACE
134
135INTERFACE
136 FUNCTION pj_inv(val, proj) bind(C,name='pj_inv')
137 IMPORT
138 TYPE(pjuv_object),VALUE :: val
139 TYPE(pj_object),VALUE :: proj
140 TYPE(pjuv_object) :: pj_inv
141 END FUNCTION pj_inv
142END INTERFACE
143
144INTERFACE
145 FUNCTION pj_geocentric_to_geodetic(a, es, point_count, point_offset, x, y, z) &
146 bind(c,name='pj_geocentric_to_geodetic')
147 IMPORT
148 REAL(kind=c_double),VALUE :: a
149 REAL(kind=c_double),VALUE :: es
150 INTEGER(kind=c_long),VALUE :: point_count
151 INTEGER(kind=c_int),VALUE :: point_offset
152 REAL(kind=c_double) :: x(*)
153 REAL(kind=c_double) :: y(*)
154 REAL(kind=c_double) :: z(*)
155 INTEGER(kind=c_int) :: pj_geocentric_to_geodetic
156 END FUNCTION pj_geocentric_to_geodetic
157END INTERFACE
158
159INTERFACE
160 FUNCTION pj_geodetic_to_geocentric(a, es, point_count, point_offset, x, y, z) &
161 bind(c,name='pj_geodetic_to_geocentric')
162 IMPORT
163 REAL(kind=c_double),VALUE :: a
164 REAL(kind=c_double),VALUE :: es
165 INTEGER(kind=c_long),VALUE :: point_count
166 INTEGER(kind=c_int),VALUE :: point_offset
167 REAL(kind=c_double) :: x(*)
168 REAL(kind=c_double) :: y(*)
169 REAL(kind=c_double) :: z(*)
170 INTEGER(kind=c_int) :: pj_geodetic_to_geocentric
171 END FUNCTION pj_geodetic_to_geocentric
172END INTERFACE
173
174INTERFACE
175 FUNCTION pj_compare_datums(srcdefn, dstdefn) bind(C,name='pj_compare_datums')
176 IMPORT
177 TYPE(pj_object),VALUE :: srcdefn
178 TYPE(pj_object),VALUE :: dstdefn
179 INTEGER(kind=c_int) :: pj_compare_datums
180 END FUNCTION pj_compare_datums
181END INTERFACE
182
183INTERFACE
184 FUNCTION pj_is_latlong(proj) bind(C,name='pj_is_latlong')
185 IMPORT
186 TYPE(pj_object),VALUE :: proj
187 INTEGER(kind=c_int) :: pj_is_latlong
188 END FUNCTION pj_is_latlong
189END INTERFACE
190
191INTERFACE
192 FUNCTION pj_is_geocent(proj) bind(C,name='pj_is_geocent')
193 IMPORT
194 TYPE(pj_object),VALUE :: proj
195 INTEGER(kind=c_int) :: pj_is_geocent
196 END FUNCTION pj_is_geocent
197END INTERFACE
198
199INTERFACE
200 FUNCTION pj_get_def(proj, options) bind(C,name='pj_get_def')
201 IMPORT
202 TYPE(pj_object),VALUE :: proj
203 INTEGER(kind=c_int) :: options
204 TYPE(c_ptr) :: pj_get_def
205 END FUNCTION pj_get_def
206END INTERFACE
207
208
209INTERFACE
210 FUNCTION pj_latlong_from_proj(proj) bind(C,name='pj_latlong_from_proj')
211 IMPORT
212 TYPE(pj_object),VALUE :: proj
213 TYPE(pj_object) :: pj_latlong_from_proj
214 END FUNCTION pj_latlong_from_proj
215END INTERFACE
216
217!INTERFACE
218! FUNCTION pj_apply_gridshift(c, i, point_count, point_offset, x, y, z) &
219! BIND(C,name='pj_apply_gridshift')
220! IMPORT
221! CHARACTER(kind=c_char) :: c(*)
222! INTEGER(kind=c_int),VALUE :: i
223! INTEGER(kind=c_long),VALUE :: point_count
224! INTEGER(kind=c_int),VALUE :: point_offset
225! REAL(kind=c_double) :: x(*)
226! REAL(kind=c_double) :: y(*)
227! REAL(kind=c_double) :: z(*)
228! INTEGER(kind=c_int) :: pj_apply_gridshift
229! END FUNCTION pj_apply_gridshift
230!END INTERFACE
231!
232!INTERFACE
233! SUBROUTINE pj_deallocate_grids() BIND(C,name='pj_deallocate_grids')
234! END SUBROUTINE pj_deallocate_grids
235!END INTERFACE
236
237INTERFACE
238 SUBROUTINE pj_free(proj) bind(C,name='pj_free')
239 IMPORT
240 TYPE(pj_object),VALUE :: proj
241 END SUBROUTINE pj_free
242END INTERFACE
243
244!void pj_pr_list(projPJ);
245!void pj_set_finder( const char *(*)(const char *) );
246!void pj_set_searchpath ( int count, const char **path );
247!projPJ pj_init(int, char **);
248!char *pj_get_def(projPJ, int);
249!void *pj_malloc(size_t);
250!void pj_dalloc(void *);
251!char *pj_strerrno(int);
252!int *pj_get_errno_ref(void);
253!const char *pj_get_release(void);
254
255CONTAINS
256
257!> Test whether the result of a pj_init_plus is a valid projection.
258!! Returns a logical value. If the second argument is provided, it
259!! checks whether they point to the same projection.
260FUNCTION pj_associated_object(arg1, arg2) RESULT(associated_)
261TYPE(pj_object),INTENT(in) :: arg1 !< projecton object to test
262TYPE(pj_object),INTENT(in),OPTIONAL :: arg2 !< optional second argument for equality test instead of validity
263LOGICAL :: associated_
264IF(PRESENT(arg2)) THEN
265 associated_ = c_associated(arg1%ptr, arg2%ptr)
266ELSE
267 associated_ = c_associated(arg1%ptr)
268ENDIF
269END FUNCTION pj_associated_object
270
271! Fortran specific version of some functions
272
273!> Fortran version of \a pj_transform proj API function.
274!! This is the Fortran version of \a pj_transform function, the array
275!! arguments are assumed-shape Fortran arrays of equal length so no
276!! array size nor offset need to be passed, see the original C API
277!! documentation for the use of the function.
278FUNCTION pj_transform_f(src, dst, x, y, z)
279TYPE(pj_object),VALUE :: src !< source coordinate system
280TYPE(pj_object),VALUE :: dst !< destination coordinate system
281REAL(kind=c_double), TARGET :: x(:) !< array of x coordinates
282REAL(kind=c_double), TARGET :: y(:) !< array of y coordinates
283REAL(kind=c_double), TARGET, OPTIONAL :: z(:) !< optional array of z coordinates
284INTEGER(kind=c_int) :: pj_transform_f
285
286REAL(kind=c_double),POINTER :: px, py, pz
287
288! a fortran pointer is required to avoid compilation errors with some
289! versions of gfortran due to x,y,z being C-incompatible assumed-shape
290! arrays
291px => x(1)
292py => y(1)
293
294IF (PRESENT(z)) THEN
295 pz => z(1)
296 pj_transform_f = pj_transform(src, dst, &
297 int(min(SIZE(x),SIZE(y),SIZE(z)), kind=c_long), 1_c_int, c_loc(px), c_loc(py), c_loc(pz))
298ELSE
299 pj_transform_f = pj_transform(src, dst, &
300 int(min(SIZE(x),SIZE(y)), kind=c_long), 1_c_int, c_loc(px), c_loc(py), c_null_ptr)
301ENDIF
302
303END FUNCTION pj_transform_f
304
305
306!> Fortran version of \a pj_datum_transform proj API function.
307!! This is the Fortran version of \a pj_datum_transform function, the array
308!! arguments are assumed-shape Fortran arrays of equal length so no
309!! array size nor offset need to be passed, see the original C API
310!! documentation for the use of the function.
311FUNCTION pj_datum_transform_f(src, dst, x, y, z)
312TYPE(pj_object),VALUE :: src !< source coordinate system
313TYPE(pj_object),VALUE :: dst !< destination coordinate system
314REAL(kind=c_double) :: x(:) !< array of x coordinates
315REAL(kind=c_double) :: y(:) !< array of y coordinates
316REAL(kind=c_double),OPTIONAL :: z(:) !< optional array of z coordinates
317INTEGER(kind=c_int) :: pj_datum_transform_f
318
319REAL(kind=c_double),POINTER :: dummyz(:)
320
321IF (PRESENT(z)) THEN
322 pj_datum_transform_f = pj_datum_transform(src, dst, &
323 int(min(SIZE(x),SIZE(y),SIZE(z)), kind=c_long), 1_c_int, x, y, z)
324ELSE
325 NULLIFY(dummyz)
326 pj_datum_transform_f = pj_datum_transform(src, dst, &
327 int(min(SIZE(x),SIZE(y)), kind=c_long), 1_c_int, x, y, dummyz)
328ENDIF
329
330END FUNCTION pj_datum_transform_f
331
332
333!> Fortran version of \a pj_geocentric_to_geodetic proj API function.
334!! This is the Fortran version of \a pj_geocentric_to_geodetic function, the array
335!! arguments are assumed-shape Fortran arrays of equal length so no
336!! array size nor offset need to be passed, see the original C API
337!! documentation for the use of the function.
338FUNCTION pj_geocentric_to_geodetic_f(a, es, x, y, z)
339REAL(kind=c_double),VALUE :: a !< Earth semi-major axis
340REAL(kind=c_double),VALUE :: es !< Earth flattening
341REAL(kind=c_double) :: x(:) !< array of x coordinates
342REAL(kind=c_double) :: y(:) !< array of y coordinates
343REAL(kind=c_double) :: z(:) !< array of z coordinates
344INTEGER(kind=c_int) :: pj_geocentric_to_geodetic_f
345
346pj_geocentric_to_geodetic_f = &
347 pj_geocentric_to_geodetic(a, es, &
348 int(min(SIZE(x),SIZE(y),SIZE(x)), kind=c_long), 1_c_int, x, y, z)
349
350END FUNCTION pj_geocentric_to_geodetic_f
351
352
353!> Fortran version of \a pj_geodetic_to_geocentric proj API function.
354!! This is the Fortran version of \a pj_geodetic_to_geocentric function, the array
355!! arguments are assumed-shape Fortran arrays of equal length so no
356!! array size nor offset need to be passed, see the original C API
357!! documentation for the use of the function.
358FUNCTION pj_geodetic_to_geocentric_f(a, es, x, y, z)
359REAL(kind=c_double),VALUE :: a !< Earth semi-major axis
360REAL(kind=c_double),VALUE :: es !< Earth flattening
361REAL(kind=c_double) :: x(:) !< array of x coordinates
362REAL(kind=c_double) :: y(:) !< array of y coordinates
363REAL(kind=c_double) :: z(:) !< array of z coordinates
364INTEGER(kind=c_int) :: pj_geodetic_to_geocentric_f
365
366pj_geodetic_to_geocentric_f = &
367 pj_geodetic_to_geocentric(a, es, &
368 int(min(SIZE(x),SIZE(y),SIZE(x)), kind=c_long), 1_c_int, x, y, z)
369
370END FUNCTION pj_geodetic_to_geocentric_f
371
372
373END MODULE proj
Fortran 2003 interface to the proj.4 https://github.com/OSGeo/proj.4 library.
Definition proj.F90:62
Object describing a cartographic projection.
Definition proj.F90:69