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
51MODULE proj
52use,INTRINSIC :: iso_c_binding
53IMPLICIT NONE
54
58TYPE,BIND(C) :: pj_object
59 PRIVATE
60 TYPE(c_ptr) :: ptr = c_null_ptr
61END TYPE pj_object
66TYPE(pj_object),PARAMETER :: pj_object_null=pj_object(c_null_ptr)
67
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
77REAL(kind=c_double),PARAMETER :: pj_rad_to_deg=57.29577951308232d0
78
82INTERFACE pj_associated
83 MODULE PROCEDURE pj_associated_object
84END INTERFACE pj_associated
85
88INTERFACE
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
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
105 TYPE(c_ptr), VALUE :: y
106 TYPE(c_ptr), VALUE :: 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
260FUNCTION pj_associated_object(arg1, arg2) RESULT(associated_)
261TYPE(pj_object),INTENT(in) :: arg1
262TYPE(pj_object),INTENT(in),OPTIONAL :: arg2
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
278FUNCTION pj_transform_f(src, dst, x, y, z)
279TYPE(pj_object),VALUE :: src
280TYPE(pj_object),VALUE :: dst
281REAL(kind=c_double), TARGET :: x(:)
282REAL(kind=c_double), TARGET :: y(:)
283REAL(kind=c_double), TARGET, OPTIONAL :: z(:)
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
311FUNCTION pj_datum_transform_f(src, dst, x, y, z)
312TYPE(pj_object),VALUE :: src
313TYPE(pj_object),VALUE :: dst
314REAL(kind=c_double) :: x(:)
315REAL(kind=c_double) :: y(:)
316REAL(kind=c_double),OPTIONAL :: z(:)
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
338FUNCTION pj_geocentric_to_geodetic_f(a, es, x, y, z)
339REAL(kind=c_double),VALUE :: a
340REAL(kind=c_double),VALUE :: es
341REAL(kind=c_double) :: x(:)
342REAL(kind=c_double) :: y(:)
343REAL(kind=c_double) :: z(:)
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
358FUNCTION pj_geodetic_to_geocentric_f(a, es, x, y, z)
359REAL(kind=c_double),VALUE :: a
360REAL(kind=c_double),VALUE :: es
361REAL(kind=c_double) :: x(:)
362REAL(kind=c_double) :: y(:)
363REAL(kind=c_double) :: z(:)
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