FortranGIS Version 3.0
Data Types | Variables
proj6 Module Reference

Fortran 2003 interface to the proj https://proj.org/ library, API version 6. More...

Data Types

type  pj_object
 Object describing a projection or a transformation. More...
 

Variables

type(pj_object), parameter pj_object_null =pj_object(C_NULL_PTR)
 Object representing a null projection or transformation.
 
type(pj_context_object), parameter pj_default_ctx =pj_context_object(C_NULL_PTR)
 Object representing the default context.
 

Detailed Description

Fortran 2003 interface to the proj https://proj.org/ library, API version 6.

The following functions and subroutines are directly interfaced to their corresponding C version, so they are undocumented here, please refer to the original proj C API documentation, e.g. at the address https://proj.org/development/reference/index.html , for their use:

Some of these functions have also a more Fortran-friendly interface explicitely documented here, with an _f appended to the name.

For an example of application of the proj6 module, please refer to the following test program, which performs a forward and backward transformation:

PROGRAM proj_test6
use,INTRINSIC :: iso_c_binding
USE proj6
IMPLICIT none
! declare projection objects
TYPE(pj_object) :: pj, pjf, pjt
! declare coordinate objects
TYPE(pj_coord_object) :: coordg, coordp, coordgc, &
vcoordg(4), vcoordp(4), vcoordgc(4)
! declare area object
TYPE(pj_area_object) :: area
INTEGER :: res
CHARACTER(len=512) :: proj_stringf, proj_stringt
! from EPSG:4326 ("latlon") to UTM
proj_stringf = 'EPSG:4326'
proj_stringt = '+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +type=crs' ! EPSG:32632
print*,'from: ',trim(proj_stringf)
print*,'to: ',trim(proj_stringt)
! initialize a projection object from two projection strings, remember //char(0)!
print*,'Defining a projection object'
pj = proj_create_crs_to_crs(pj_default_ctx, trim(proj_stringf)//char(0), &
trim(proj_stringt)//char(0), area)
print*,'Converting through a pj_coord object'
! define a coordinate set, %z and %t keep the default value of 0;
! EPSG:4326 expects values in degrees, lat lon order
coordg%x = 45.0d0
coordg%y = 11.0d0
! make a forward projection to UTM, result in meters x,y order
coordp = proj_trans(pj, pj_fwd, coordg)
! return back to polar coordinates in degrees
coordgc = proj_trans(pj, pj_inv, coordp)
print*,'Original:',coordg
print*,'Projected:',coordp
print*,'Returned:',coordgc
! check whether they are the same within 10-8 degrees
IF (abs(coordgc%x-coordg%x) > 1.0d-8 .OR. abs(coordgc%y-coordg%y) > 1.0d-8) THEN
print*,'Error pj_inv*pj_fwd failed or /= identity'
CALL exit(1)
ENDIF
pj = proj_destroy(pj)
! initialize a projection transformation object from two projection objects
pjf = proj_create(pj_default_ctx, trim(proj_stringf)//char(0))
IF (.NOT.proj_associated(pjf)) THEN
print*,'Failing to define projection object from ',trim(proj_stringf)
CALL print_error(pj_default_ctx)
CALL exit(1)
ENDIF
CALL print_info(proj_pj_info(pjf))
pjt = proj_create(pj_default_ctx, trim(proj_stringt)//char(0))
IF (.NOT.proj_associated(pjt)) THEN
print*,'Failing to define projection object from ',trim(proj_stringt)
CALL print_error(pj_default_ctx)
CALL exit(1)
ENDIF
CALL print_info(proj_pj_info(pjt))
print*,'Defining a projection object'
pj = proj_create_crs_to_crs_from_pj(pj_default_ctx, pjf, pjt, area, c_null_ptr)
IF (.NOT.proj_associated(pjt)) THEN
print*,'Failing to define projection object'
CALL print_error(pj_default_ctx)
CALL exit(1)
ENDIF
print*,'Converting through a pj_coord array object'
! define a coordinate set, values in degree, %z and %t keep the
! default value of 0
vcoordg(:)%x = (/45.0d0,45.0d0,46.0d0,46.0d0/)
vcoordg(:)%y = (/11.0d0,12.0d0,11.0d0,12.0d0/)
vcoordp = vcoordg ! array transformation overwrites
! make a forward projection to UTM, result in meters
res = proj_trans_f(pj, pj_fwd, vcoordp)
IF (res /= 0) THEN
print*,'Failure in forward array transformation',res
CALL print_error(pj_default_ctx)
ENDIF
vcoordgc = vcoordp ! array transformation overwrites
! return back to polar coordinates in degrees
res = proj_trans_f(pj, pj_inv, vcoordgc)
IF (res /= 0) THEN
print*,'Failure in inverse array transformation',res
CALL print_error(pj_default_ctx)
CALL exit(1)
ENDIF
print*,'Original:',vcoordg%x,vcoordg%y
print*,'Projected:',vcoordp%x,vcoordp%y
print*,'Returned:',vcoordgc%x,vcoordgc%y
! check whether they are the same within 10-8 degrees
IF (maxval(abs(vcoordgc%x-vcoordg%x)) > 1.0d-8 .OR. &
maxval(abs(vcoordgc%y-vcoordg%y)) > 1.0d-8) THEN
print*,'Error pj_inv*pj_fwd failed or /= identity'
CALL exit(1)
ENDIF
pj = proj_destroy(pj)
CONTAINS
SUBROUTINE print_info(info)
TYPE(pj_proj_info_object),INTENT(in) :: info
print*,'==== Object information ===='
print*,'id:',trim(strtofchar(info%id,256))
print*,'desc:',trim(strtofchar(info%description,256))
print*,'def:',trim(strtofchar(info%definition,256))
print*,'has inv:',info%has_inverse
print*,'accuracy:',info%accuracy
print*,'============================'
END SUBROUTINE print_info
SUBROUTINE print_error(ctx)
TYPE(pj_context_object),INTENT(in) :: ctx
IF (proj_context_errno(ctx) /= 0) THEN
! switch to proj_context_errno_string
print*,trim(strtofchar(proj_errno_string(proj_context_errno(ctx)), 256))
ENDIF
END SUBROUTINE print_error
END 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