FortranGIS  Version2.5
proj_test.F90
1 PROGRAM proj_test
2 use,INTRINSIC :: iso_c_binding
3 USE fortranc
4 USE proj
5 IMPLICIT none
6 
7 ! declare projection objects
8 TYPE(pj_object) :: pj, pjll
9 ! declare coordinate objects
10 TYPE(pjuv_object) :: coordg, coordp, coordgc
11 REAL(kind=c_double) :: x(4), y(4), z(4), xorig(4), yorig(4)
12 INTEGER :: res
13 CHARACTER(len=512) :: proj_string
14 
15 proj_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)!
18 print*,trim(proj_string)
19 pj = pj_init_plus(trim(proj_string)//char(0))
20 print*,'Associated:',pj_associated(pj)
21 IF (.NOT.pj_associated(pj)) CALL exit(1)
22 print*,'Geographic:',pj_is_latlong(pj)
23 
24 ! build a corresponding latlon projection
25 print*,'Corresponding latlon projection'
26 pjll = pj_latlong_from_proj(pj)
27 IF (.NOT.pj_associated(pjll)) CALL exit(1)
28 print*,'Associated:',pj_associated(pjll)
29 print*,'Geographic:',pj_is_latlong(pjll)
30 ! retrieve the projection string defined by proj, the trick described
31 ! in the fortranc library is used
32 proj_string = strtofchar(pj_get_def(pjll, 0), len(proj_string))
33 print*,trim(proj_string)
34 
35 print*,'Converting through a pjuv object'
36 ! define a coordinate set, it is latlong so convert to radians
37 coordg = pjuv_object(11.0d0*pj_deg_to_rad, 45.0d0*pj_deg_to_rad)
38 
39 ! make a forward projection, result in meters
40 coordp = pj_fwd(coordg, pj)
41 ! return back to polar coordinates in radians
42 coordgc = pj_inv(coordp, pj)
43 print*,'Original:',coordg
44 print*,'Projected:',coordp
45 print*,'Returned:',coordgc
46 
47 ! check whether they are the same within 10-6 radians
48 IF (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)
51 ENDIF
52 
53 xorig(:) = (/11.0d0,12.0d0,11.0d0,12.0d0/)*pj_deg_to_rad
54 yorig(:) = (/45.0d0,45.0d0,46.0d0,46.0d0/)*pj_deg_to_rad
55 
56 print*,'Converting through pj_transform_f with z'
57 x(:) = xorig(:)
58 y(:) = yorig(:)
59 z(:) = 0.0d0
60 print*,'Original:',x(1),y(1)
61 res = pj_transform_f(pjll, pj, x, y, z)
62 print*,'Projected:',x(1),y(1)
63 IF (res == 0) res = pj_transform_f(pj, pjll, x, y, z)
64 IF (res == 0) THEN
65  print*,'Returned:',x(1),y(1)
66 ELSE
67  print*,'pj_transform with z failed'
68  CALL exit(1)
69 ENDIF
70 IF (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)
73 ENDIF
74 
75 print*,'Converting through pj_transform_f without z'
76 x(:) = xorig(:)
77 y(:) = yorig(:)
78 print*,'Original:',x(1),y(1)
79 res = pj_transform_f(pjll, pj, x, y)
80 print*,'Projected:',x(1),y(1)
81 IF (res == 0) res = pj_transform_f(pj, pjll, x, y)
82 IF (res == 0) THEN
83  print*,'Returned:',x(1),y(1)
84 ELSE
85  print*,'pj_transform without z failed'
86  CALL exit(1)
87 ENDIF
88 IF (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)
91 ENDIF
92 
93 CALL pj_free(pj)
94 CALL pj_free(pjll)
95 
96 !x(1) = coordg%u
97 !y(1) = coordg%v
98 !z(1) = 0.0D0
99 !res = pj_transform(pjll, pj, x, y, z)
100 !PRINT*,x,y,res
101 
102 END PROGRAM proj_test
103 
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 coordinate pair.
Definition: proj.F90:82
Object describing a cartographic projection.
Definition: proj.F90:69
Utility module for supporting Fortran 2003 C language interface module.
Definition: fortranc.F90:100
Convert a null-terminated C string into a Fortran CHARACTER variable of the proper length...
Definition: fortranc.F90:165
Test whether an opaque object is valid.
Definition: proj.F90:93