proj_test.F90

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