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 
51 MODULE proj
52 use,INTRINSIC :: iso_c_binding
53 IMPLICIT NONE
54 
58 TYPE,BIND(C) :: pj_object
59  PRIVATE
60  TYPE(c_ptr) :: ptr = c_null_ptr
61 END TYPE pj_object
62 
66 TYPE(pj_object),PARAMETER :: pj_object_null=pj_object(c_null_ptr)
67 
71 TYPE,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)
74 END TYPE pjuv_object
75 
76 REAL(kind=c_double),PARAMETER :: pj_deg_to_rad=.0174532925199432958d0
77 REAL(kind=c_double),PARAMETER :: pj_rad_to_deg=57.29577951308232d0
78 
82 INTERFACE pj_associated
83  MODULE PROCEDURE pj_associated_object
84 END INTERFACE pj_associated
85 
88 INTERFACE
89  FUNCTION pj_init_plus(name) bind(C,name='pj_init_plus')
90  IMPORT
91  CHARACTER(kind=c_char) :: name(*)
92  TYPE(pj_object) :: pj_init_plus
93  END FUNCTION pj_init_plus
94 END INTERFACE
95 
96 INTERFACE
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
105  TYPE(c_ptr), VALUE :: y
106  TYPE(c_ptr), VALUE :: z
107  INTEGER(kind=c_int) :: pj_transform
108  END FUNCTION pj_transform
109 END INTERFACE
110 
111 INTERFACE
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
124 END INTERFACE
125 
126 INTERFACE
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
133 END INTERFACE
134 
135 INTERFACE
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
142 END INTERFACE
143 
144 INTERFACE
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
157 END INTERFACE
158 
159 INTERFACE
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
172 END INTERFACE
173 
174 INTERFACE
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
181 END INTERFACE
182 
183 INTERFACE
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
189 END INTERFACE
190 
191 INTERFACE
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
197 END INTERFACE
198 
199 INTERFACE
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
206 END INTERFACE
207 
208 
209 INTERFACE
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
215 END 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 
237 INTERFACE
238  SUBROUTINE pj_free(proj) bind(C,name='pj_free')
239  IMPORT
240  TYPE(pj_object),VALUE :: proj
241  END SUBROUTINE pj_free
242 END 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 
255 CONTAINS
256 
260 FUNCTION pj_associated_object(arg1, arg2) RESULT(associated_)
261 TYPE(pj_object),INTENT(in) :: arg1
262 TYPE(pj_object),INTENT(in),OPTIONAL :: arg2
263 LOGICAL :: associated_
264 IF(PRESENT(arg2)) THEN
265  associated_ = c_associated(arg1%ptr, arg2%ptr)
266 ELSE
267  associated_ = c_associated(arg1%ptr)
268 ENDIF
269 END FUNCTION pj_associated_object
270 
271 ! Fortran specific version of some functions
272 
278 FUNCTION pj_transform_f(src, dst, x, y, z)
279 TYPE(pj_object),VALUE :: src
280 TYPE(pj_object),VALUE :: dst
281 REAL(kind=c_double), TARGET :: x(:)
282 REAL(kind=c_double), TARGET :: y(:)
283 REAL(kind=c_double), TARGET, OPTIONAL :: z(:)
284 INTEGER(kind=c_int) :: pj_transform_f
285 
286 REAL(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
291 px => x(1)
292 py => y(1)
293 
294 IF (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))
298 ELSE
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)
301 ENDIF
302 
303 END FUNCTION pj_transform_f
304 
305 
311 FUNCTION pj_datum_transform_f(src, dst, x, y, z)
312 TYPE(pj_object),VALUE :: src
313 TYPE(pj_object),VALUE :: dst
314 REAL(kind=c_double) :: x(:)
315 REAL(kind=c_double) :: y(:)
316 REAL(kind=c_double),OPTIONAL :: z(:)
317 INTEGER(kind=c_int) :: pj_datum_transform_f
318 
319 REAL(kind=c_double),POINTER :: dummyz(:)
320 
321 IF (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)
324 ELSE
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)
328 ENDIF
329 
330 END FUNCTION pj_datum_transform_f
331 
332 
338 FUNCTION pj_geocentric_to_geodetic_f(a, es, x, y, z)
339 REAL(kind=c_double),VALUE :: a
340 REAL(kind=c_double),VALUE :: es
341 REAL(kind=c_double) :: x(:)
342 REAL(kind=c_double) :: y(:)
343 REAL(kind=c_double) :: z(:)
344 INTEGER(kind=c_int) :: pj_geocentric_to_geodetic_f
345 
346 pj_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 
350 END FUNCTION pj_geocentric_to_geodetic_f
351 
352 
358 FUNCTION pj_geodetic_to_geocentric_f(a, es, x, y, z)
359 REAL(kind=c_double),VALUE :: a
360 REAL(kind=c_double),VALUE :: es
361 REAL(kind=c_double) :: x(:)
362 REAL(kind=c_double) :: y(:)
363 REAL(kind=c_double) :: z(:)
364 INTEGER(kind=c_int) :: pj_geodetic_to_geocentric_f
365 
366 pj_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 
370 END FUNCTION pj_geodetic_to_geocentric_f
371 
372 
373 END MODULE proj
Test whether an opaque object is valid.
Definition: proj.F90:93
Initialize a projection from a string.
Definition: proj.F90:100
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
Object describing a coordinate pair.
Definition: proj.F90:82