2use,
INTRINSIC :: iso_c_binding
8TYPE(pj_object) :: pj, pjll
10TYPE(pjuv_object) :: coordg, coordp, coordgc
11REAL(kind=c_double) :: x(4), y(4), z(4), xorig(4), yorig(4)
13CHARACTER(len=512) :: proj_string
16 '+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0'
18print*,trim(proj_string)
19pj = pj_init_plus(trim(proj_string)//char(0))
20print*,
'Associated:',pj_associated(pj)
21IF (.NOT.pj_associated(pj))
CALL exit(1)
22print*,
'Geographic:',pj_is_latlong(pj)
25print*,
'Corresponding latlon projection'
26pjll = pj_latlong_from_proj(pj)
27IF (.NOT.pj_associated(pjll))
CALL exit(1)
28print*,
'Associated:',pj_associated(pjll)
29print*,
'Geographic:',pj_is_latlong(pjll)
32proj_string =
strtofchar(pj_get_def(pjll, 0), len(proj_string))
33print*,trim(proj_string)
35print*,
'Converting through a pjuv object'
37coordg = pjuv_object(11.0d0*pj_deg_to_rad, 45.0d0*pj_deg_to_rad)
40coordp = pj_fwd(coordg, pj)
42coordgc = pj_inv(coordp, pj)
43print*,
'Original:',coordg
44print*,
'Projected:',coordp
45print*,
'Returned:',coordgc
48IF (abs(coordgc%u-coordg%u) > 1.0d-6 .OR. abs(coordgc%v-coordg%v) > 1.0d-6)
THEN
49 print*,
'Error pj_inv*pj_fwd failed or /= identity'
53xorig(:) = (/11.0d0,12.0d0,11.0d0,12.0d0/)*pj_deg_to_rad
54yorig(:) = (/45.0d0,45.0d0,46.0d0,46.0d0/)*pj_deg_to_rad
56print*,
'Converting through pj_transform_f with z'
60print*,
'Original:',x(1),y(1)
61res = pj_transform_f(pjll, pj, x, y, z)
62print*,
'Projected:',x(1),y(1)
63IF (res == 0) res = pj_transform_f(pj, pjll, x, y, z)
65 print*,
'Returned:',x(1),y(1)
67 print*,
'pj_transform with z failed'
70IF (maxval(abs(xorig-x)) > 1.0d-6 .OR. maxval(abs(yorig-y)) > 1.0d-6)
THEN
71 print*,
'Error pj_transform*pj_transform**-1 /= identity'
75print*,
'Converting through pj_transform_f without z'
78print*,
'Original:',x(1),y(1)
79res = pj_transform_f(pjll, pj, x, y)
80print*,
'Projected:',x(1),y(1)
81IF (res == 0) res = pj_transform_f(pj, pjll, x, y)
83 print*,
'Returned:',x(1),y(1)
85 print*,
'pj_transform without z failed'
88IF (maxval(abs(xorig-x)) > 1.0d-6 .OR. maxval(abs(yorig-y)) > 1.0d-6)
THEN
89 print*,
'Error pj_transform*pj_transform**-1 /= identity'
Convert a null-terminated C string into a Fortran CHARACTER variable of the proper length.
Utility module for supporting Fortran 2003 C language interface module.
Fortran 2003 interface to the proj.4 https://github.com/OSGeo/proj.4 library.