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