2 use,
INTRINSIC :: iso_c_binding
8 TYPE(pj_object) :: pj, pjf, pjt
10 TYPE(pj_coord_object) :: coordg, coordp, coordgc, &
11 vcoordg(4), vcoordp(4), vcoordgc(4)
13 TYPE(pj_area_object) :: area
16 CHARACTER(len=512) :: proj_stringf, proj_stringt
19 proj_stringf =
'EPSG:4326'
20 proj_stringt =
'+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +type=crs'
21 print*,
'from: ',trim(proj_stringf)
22 print*,
'to: ',trim(proj_stringt)
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)
29 print*,
'Converting through a pj_coord object'
36 coordp = proj_trans(pj, pj_fwd, coordg)
38 coordgc = proj_trans(pj, pj_inv, coordp)
40 print*,
'Original:',coordg
41 print*,
'Projected:',coordp
42 print*,
'Returned:',coordgc
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'
53 pjf = proj_create(pj_default_ctx, trim(proj_stringf)//char(0))
55 print*,
'Failing to define projection object from ',trim(proj_stringf)
56 CALL print_error(pj_default_ctx)
59 CALL print_info(proj_pj_info(pjf))
61 pjt = proj_create(pj_default_ctx, trim(proj_stringt)//char(0))
63 print*,
'Failing to define projection object from ',trim(proj_stringt)
64 CALL print_error(pj_default_ctx)
67 CALL print_info(proj_pj_info(pjt))
69 print*,
'Defining a projection object'
70 pj = proj_create_crs_to_crs_from_pj(pj_default_ctx, pjf, pjt, area, c_null_ptr)
72 print*,
'Failing to define projection object'
73 CALL print_error(pj_default_ctx)
77 print*,
'Converting through a pj_coord array object'
80 vcoordg(:)%x = (/45.0d0,45.0d0,46.0d0,46.0d0/)
81 vcoordg(:)%y = (/11.0d0,12.0d0,11.0d0,12.0d0/)
85 res = proj_trans_f(pj, pj_fwd, vcoordp)
87 print*,
'Failure in forward array transformation',res
88 CALL print_error(pj_default_ctx)
93 res = proj_trans_f(pj, pj_inv, vcoordgc)
95 print*,
'Failure in inverse array transformation',res
96 CALL print_error(pj_default_ctx)
100 print*,
'Original:',vcoordg%x,vcoordg%y
101 print*,
'Projected:',vcoordp%x,vcoordp%y
102 print*,
'Returned:',vcoordgc%x,vcoordgc%y
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'
111 pj = proj_destroy(pj)
115 SUBROUTINE print_info(info)
116 TYPE(pj_proj_info_object),
INTENT(in) :: info
118 print*,
'==== Object information ===='
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*,
'============================'
126 END SUBROUTINE print_info
128 SUBROUTINE print_error(ctx)
129 TYPE(pj_context_object),
INTENT(in) :: ctx
131 IF (proj_context_errno(ctx) /= 0)
THEN
133 print*,trim(
strtofchar(proj_errno_string(proj_context_errno(ctx)), 256))
136 END SUBROUTINE print_error
138 END PROGRAM proj_test6
Convert a null-terminated C string into a Fortran CHARACTER variable of the proper length.
Test whether an opaque object is valid.
Utility module for supporting Fortran 2003 C language interface module.
Fortran 2003 interface to the proj https://proj.org/ library, API version 6.