FortranGIS Version 3.0
proj6_test.F90
1PROGRAM proj_test6
2use,INTRINSIC :: iso_c_binding
3USE fortranc
4USE proj6
5IMPLICIT none
6
7! declare projection objects
8TYPE(pj_object) :: pj, pjf, pjt
9! declare coordinate objects
10TYPE(pj_coord_object) :: coordg, coordp, coordgc, &
11 vcoordg(4), vcoordp(4), vcoordgc(4)
12! declare area object
13TYPE(pj_area_object) :: area
14
15INTEGER :: res
16CHARACTER(len=512) :: proj_stringf, proj_stringt
17
18! from EPSG:4326 ("latlon") to UTM
19proj_stringf = 'EPSG:4326'
20proj_stringt = '+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +type=crs' ! EPSG:32632
21print*,'from: ',trim(proj_stringf)
22print*,'to: ',trim(proj_stringt)
23
24! initialize a projection object from two projection strings, remember //char(0)!
25print*,'Defining a projection object'
26pj = proj_create_crs_to_crs(pj_default_ctx, trim(proj_stringf)//char(0), &
27 trim(proj_stringt)//char(0), area)
28
29print*,'Converting through a pj_coord object'
30! define a coordinate set, %z and %t keep the default value of 0;
31! EPSG:4326 expects values in degrees, lat lon order
32coordg%x = 45.0d0
33coordg%y = 11.0d0
34
35! make a forward projection to UTM, result in meters x,y order
36coordp = proj_trans(pj, pj_fwd, coordg)
37! return back to polar coordinates in degrees
38coordgc = proj_trans(pj, pj_inv, coordp)
39
40print*,'Original:',coordg
41print*,'Projected:',coordp
42print*,'Returned:',coordgc
43
44! check whether they are the same within 10-8 degrees
45IF (abs(coordgc%x-coordg%x) > 1.0d-8 .OR. abs(coordgc%y-coordg%y) > 1.0d-8) THEN
46 print*,'Error pj_inv*pj_fwd failed or /= identity'
47 CALL exit(1)
48ENDIF
49
50pj = proj_destroy(pj)
51
52! initialize a projection transformation object from two projection objects
53pjf = proj_create(pj_default_ctx, trim(proj_stringf)//char(0))
54IF (.NOT.proj_associated(pjf)) THEN
55 print*,'Failing to define projection object from ',trim(proj_stringf)
56 CALL print_error(pj_default_ctx)
57 CALL exit(1)
58ENDIF
59CALL print_info(proj_pj_info(pjf))
60
61pjt = proj_create(pj_default_ctx, trim(proj_stringt)//char(0))
62IF (.NOT.proj_associated(pjt)) THEN
63 print*,'Failing to define projection object from ',trim(proj_stringt)
64 CALL print_error(pj_default_ctx)
65 CALL exit(1)
66ENDIF
67CALL print_info(proj_pj_info(pjt))
68
69print*,'Defining a projection object'
70pj = proj_create_crs_to_crs_from_pj(pj_default_ctx, pjf, pjt, area, c_null_ptr)
71IF (.NOT.proj_associated(pjt)) THEN
72 print*,'Failing to define projection object'
73 CALL print_error(pj_default_ctx)
74 CALL exit(1)
75ENDIF
76
77print*,'Converting through a pj_coord array object'
78! define a coordinate set, values in degree, %z and %t keep the
79! default value of 0
80vcoordg(:)%x = (/45.0d0,45.0d0,46.0d0,46.0d0/)
81vcoordg(:)%y = (/11.0d0,12.0d0,11.0d0,12.0d0/)
82
83vcoordp = vcoordg ! array transformation overwrites
84! make a forward projection to UTM, result in meters
85res = proj_trans_f(pj, pj_fwd, vcoordp)
86IF (res /= 0) THEN
87 print*,'Failure in forward array transformation',res
88 CALL print_error(pj_default_ctx)
89ENDIF
90vcoordgc = vcoordp ! array transformation overwrites
91
92! return back to polar coordinates in degrees
93res = proj_trans_f(pj, pj_inv, vcoordgc)
94IF (res /= 0) THEN
95 print*,'Failure in inverse array transformation',res
96 CALL print_error(pj_default_ctx)
97 CALL exit(1)
98ENDIF
99
100print*,'Original:',vcoordg%x,vcoordg%y
101print*,'Projected:',vcoordp%x,vcoordp%y
102print*,'Returned:',vcoordgc%x,vcoordgc%y
103
104! check whether they are the same within 10-8 degrees
105IF (maxval(abs(vcoordgc%x-vcoordg%x)) > 1.0d-8 .OR. &
106 maxval(abs(vcoordgc%y-vcoordg%y)) > 1.0d-8) THEN
107 print*,'Error pj_inv*pj_fwd failed or /= identity'
108 CALL exit(1)
109ENDIF
110
111pj = proj_destroy(pj)
112
113CONTAINS
114
115SUBROUTINE print_info(info)
116TYPE(pj_proj_info_object),INTENT(in) :: info
117
118print*,'==== Object information ===='
119print*,'id:',trim(strtofchar(info%id,256))
120print*,'desc:',trim(strtofchar(info%description,256))
121print*,'def:',trim(strtofchar(info%definition,256))
122print*,'has inv:',info%has_inverse
123print*,'accuracy:',info%accuracy
124print*,'============================'
125
126END SUBROUTINE print_info
127
128SUBROUTINE print_error(ctx)
129TYPE(pj_context_object),INTENT(in) :: ctx
130
131IF (proj_context_errno(ctx) /= 0) THEN
132 ! switch to proj_context_errno_string
133 print*,trim(strtofchar(proj_errno_string(proj_context_errno(ctx)), 256))
134ENDIF
135
136END SUBROUTINE print_error
137
138END PROGRAM proj_test6
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 https://proj.org/ library, API version 6.
Definition: proj6.F90:63