proj Module Reference

Fortran 2003 interface to the proj.4 https://github.com/OSGeo/proj.4 library. More...

Data Types
interface	pj_associated
	Test whether an opaque object is valid. More...
interface	pj_init_plus
	Initialize a projection from a string. More...
type	pj_object
	Object describing a cartographic projection. More...
type	pjuv_object
	Object describing a coordinate pair. More...

Functions/Subroutines
logical function	pj_associated_object (arg1, arg2)
	Test whether the result of a pj_init_plus is a valid projection. More...
integer(kind=c_int) function	pj_transform_f (src, dst, x, y, z)
	Fortran version of pj_transform proj API function. More...
integer(kind=c_int) function	pj_datum_transform_f (src, dst, x, y, z)
	Fortran version of pj_datum_transform proj API function. More...
integer(kind=c_int) function	pj_geocentric_to_geodetic_f (a, es, x, y, z)
	Fortran version of pj_geocentric_to_geodetic proj API function. More...
integer(kind=c_int) function	pj_geodetic_to_geocentric_f (a, es, x, y, z)
	Fortran version of pj_geodetic_to_geocentric proj API function. More...

Variables
type(pj_object), parameter	pj_object_null =pj_object(C_NULL_PTR)
	Object representing a null cartographic projection. More...
real(kind=c_double), parameter	pj_deg_to_rad =.0174532925199432958D0
	equivalent to the C api symbol DEG_TO_RAD More...
real(kind=c_double), parameter	pj_rad_to_deg =57.29577951308232D0
	equivalent to the C api symbol RAD_TO_DEG More...

Detailed Description

Fortran 2003 interface to the proj.4 https://github.com/OSGeo/proj.4 library.

The following functions are directly interfaced to their corresponding C version, so they are undocumented here, please refer to the original gdal C API documentation, e.g. at the address https://github.com/OSGeo/proj.4/wiki/ProjAPI , for their use:

• pj_init_plus() -> FUNCTION pj_init_plus()
• pj_transform()() -> FUNCTION pj_transform()()
• pj_datum_transform() -> FUNCTION pj_datum_transform()
• pj_fwd() -> FUNCTION pj_fwd()
• pj_inv() -> FUNCTION pj_inv()
• pj_geocentric_to_geodetic() -> FUNCTION pj_geocentric_to_geodetic()
• pj_geodetic_to_geocentric() -> FUNCTION pj_geodetic_to_geocentric()
• pj_compare_datums() -> FUNCTION pj_compare_datums()
• pj_is_latlong() -> FUNCTION pj_is_latlong()
• pj_is_geocent() -> FUNCTION pj_is_geocent()
• pj_get_def() -> FUNCTION pj_get_def()
• pj_latlong_from_proj() -> FUNCTION pj_latlong_from_proj()
• pj_free() -> SUBROUTINE pj_free()

Notice that, if relevant, the result of functions returning an integer has to be interpreted as 0=false, nonzero=true or 0=ok, nonzero=error.

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 proj module, please refer to the following test program, which performs a forward and backward transformation:

PROGRAM proj_test
use,INTRINSIC :: iso_c_binding
USE fortranc
USE proj
IMPLICIT none

! declare projection objects
TYPE(pj_object) :: pj, pjll
! declare coordinate objects
TYPE(pjuv_object) :: coordg, coordp, coordgc
REAL(kind=c_double) :: x(4), y(4), z(4), xorig(4), yorig(4)
INTEGER :: res
CHARACTER(len=512) :: proj_string

proj_string = &
 '+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0'
! initialize a projection object from a string, remember //char(0)!
print*,trim(proj_string)
pj = pj_init_plus(trim(proj_string)//char(0))
print*,'Associated:',pj_associated(pj)
IF (.NOT.pj_associated(pj)) CALL exit(1)
print*,'Geographic:',pj_is_latlong(pj)

! build a corresponding latlon projection
print*,'Corresponding latlon projection'
pjll = pj_latlong_from_proj(pj)
IF (.NOT.pj_associated(pjll)) CALL exit(1)
print*,'Associated:',pj_associated(pjll)
print*,'Geographic:',pj_is_latlong(pjll)
! retrieve the projection string defined by proj, the trick described
! in the fortranc library is used
proj_string = strtofchar(pj_get_def(pjll, 0), len(proj_string))
print*,trim(proj_string)

print*,'Converting through a pjuv object'
! define a coordinate set, it is latlong so convert to radians
coordg = pjuv_object(11.0d0*pj_deg_to_rad, 45.0d0*pj_deg_to_rad)

! make a forward projection, result in meters
coordp = pj_fwd(coordg, pj)
! return back to polar coordinates in radians
coordgc = pj_inv(coordp, pj)
print*,'Original:',coordg
print*,'Projected:',coordp
print*,'Returned:',coordgc

! check whether they are the same within 10-6 radians
IF (abs(coordgc%u-coordg%u) > 1.0d-6 .OR. abs(coordgc%v-coordg%v) > 1.0d-6) THEN
  print*,'Error pj_inv*pj_fwd failed or /= identity'
  CALL exit(1)
ENDIF

xorig(:) = (/11.0d0,12.0d0,11.0d0,12.0d0/)*pj_deg_to_rad
yorig(:) = (/45.0d0,45.0d0,46.0d0,46.0d0/)*pj_deg_to_rad

print*,'Converting through pj_transform_f with z'
x(:) = xorig(:)
y(:) = yorig(:)
z(:) = 0.0d0
print*,'Original:',x(1),y(1)
res = pj_transform_f(pjll, pj, x, y, z)
print*,'Projected:',x(1),y(1)
IF (res == 0) res = pj_transform_f(pj, pjll, x, y, z)
IF (res == 0) THEN
  print*,'Returned:',x(1),y(1)
ELSE
  print*,'pj_transform with z failed'
  CALL exit(1)
ENDIF
IF (maxval(abs(xorig-x)) > 1.0d-6 .OR. maxval(abs(yorig-y)) > 1.0d-6) THEN
  print*,'Error pj_transform*pj_transform**-1 /= identity'
  CALL exit(1)
ENDIF

print*,'Converting through pj_transform_f without z'
x(:) = xorig(:)
y(:) = yorig(:)
print*,'Original:',x(1),y(1)
res = pj_transform_f(pjll, pj, x, y)
print*,'Projected:',x(1),y(1)
IF (res == 0) res = pj_transform_f(pj, pjll, x, y)
IF (res == 0) THEN
  print*,'Returned:',x(1),y(1)
ELSE
  print*,'pj_transform without z failed'
  CALL exit(1)
ENDIF
IF (maxval(abs(xorig-x)) > 1.0d-6 .OR. maxval(abs(yorig-y)) > 1.0d-6) THEN
  print*,'Error pj_transform*pj_transform**-1 /= identity'
  CALL exit(1)
ENDIF

CALL pj_free(pj)
CALL pj_free(pjll)

!x(1) = coordg%u
!y(1) = coordg%v
!z(1) = 0.0D0
!res = pj_transform(pjll, pj, x, y, z)
!PRINT*,x,y,res

END PROGRAM proj_test