FortranGIS Version 3.0
proj_test.F90
1PROGRAM proj_test
2use,INTRINSIC :: iso_c_binding
3USE fortranc
4USE proj
5IMPLICIT none
6
7! declare projection objects
8TYPE(pj_object) :: pj, pjll
9! declare coordinate objects
10TYPE(pjuv_object) :: coordg, coordp, coordgc
11REAL(kind=c_double) :: x(4), y(4), z(4), xorig(4), yorig(4)
12INTEGER :: res
13CHARACTER(len=512) :: proj_string
14
15proj_string = &
16 '+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0'
17! initialize a projection object from a string, remember //char(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)
23
24! build a corresponding latlon projection
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)
30! retrieve the projection string defined by proj, the trick described
31! in the fortranc library is used
32proj_string = strtofchar(pj_get_def(pjll, 0), len(proj_string))
33print*,trim(proj_string)
34
35print*,'Converting through a pjuv object'
36! define a coordinate set, it is latlong so convert to radians
37coordg = pjuv_object(11.0d0*pj_deg_to_rad, 45.0d0*pj_deg_to_rad)
38
39! make a forward projection, result in meters
40coordp = pj_fwd(coordg, pj)
41! return back to polar coordinates in radians
42coordgc = pj_inv(coordp, pj)
43print*,'Original:',coordg
44print*,'Projected:',coordp
45print*,'Returned:',coordgc
46
47! check whether they are the same within 10-6 radians
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'
50 CALL exit(1)
51ENDIF
52
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
55
56print*,'Converting through pj_transform_f with z'
57x(:) = xorig(:)
58y(:) = yorig(:)
59z(:) = 0.0d0
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)
64IF (res == 0) THEN
65 print*,'Returned:',x(1),y(1)
66ELSE
67 print*,'pj_transform with z failed'
68 CALL exit(1)
69ENDIF
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'
72 CALL exit(1)
73ENDIF
74
75print*,'Converting through pj_transform_f without z'
76x(:) = xorig(:)
77y(:) = yorig(:)
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)
82IF (res == 0) THEN
83 print*,'Returned:',x(1),y(1)
84ELSE
85 print*,'pj_transform without z failed'
86 CALL exit(1)
87ENDIF
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'
90 CALL exit(1)
91ENDIF
92
93CALL pj_free(pj)
94CALL pj_free(pjll)
95
96END PROGRAM proj_test
97
Convert a null-terminated C string into a Fortran CHARACTER variable of the proper length.
Definition: fortranc.F90:174
Utility module for supporting Fortran 2003 C language interface module.
Definition: fortranc.F90:103
Fortran 2003 interface to the proj.4 https://github.com/OSGeo/proj.4 library.
Definition: proj.F90:62