<FONT face="Default Sans Serif,Verdana,Arial,Helvetica,sans-serif" size=2>Hi!<br><br>As promised:<br>For those who want to call Proj from within their Fortran program, here is a solution!<br><br>In your source-directory place a copy of module libmyp.f90:<br><br>MODULE LibMyP<br>!<br>! Module for interfacing between Proj-4 library libproj.a<br>! and Fortran. So far only pj_init_plus and pj_transform are supported.<br>! More routines can be added in a similar way.<br>! Sample use:<br>!<br>! g95 libmyp.f90 testproj.f90 -L. -lproj -o testproj.exe<br>!<br>! Formulate your compile statement with care!<br>! FIRST the sources, THEN the library-path and library!!<br>!<br>USE ISO_C_BINDING<br>IMPLICIT NONE<br>PRIVATE<br>PUBLIC :: pj_init_plus, pj_transform<br>!<br>INTERFACE<br>   FUNCTION pj_init_plus(name) BIND(C,NAME='pj_init_plus')<br>      USE ISO_C_BINDING<br>      IMPLICIT NONE<br>      TYPE(C_PTR) :: pj_init_plus<br>      CHARACTER(KIND=C_CHAR) :: name(*)<br>   END FUNCTION pj_init_plus<br><br>   FUNCTION pj_transform(src, dst, point_count, point_offset, &<br>     &   x, y, z) BIND(C,NAME='pj_transform')<br>      USE ISO_C_BINDING<br>      IMPLICIT NONE<br>      TYPE(C_PTR), VALUE :: src, dst<br>      INTEGER(C_LONG), VALUE :: point_count<br>      INTEGER(C_INT), VALUE :: point_offset<br>      INTEGER(C_INT) :: pj_transform<br>      REAL(C_DOUBLE) :: x, y, z<br>   END FUNCTION pj_transform<br>END INTERFACE<br>!<br>END MODULE LibMyP<br><br><br><br>Now the following program can be tested:<br><br><br>   PROGRAM TestProj<br>!<br>   USE ISO_C_BINDING<br>   USE libmyp<br>   IMPLICIT NONE<br>!<br>   type(C_PTR) :: SystemInPtr,SystemOutPtr<br>   real(C_DOUBLE) :: x, y, z<br>   INTEGER(C_INT) :: MyResult<br>!<br>! Initialize the coordinate systems<br>!<br>   SystemInPtr = pj_init_plus('+proj=longlat +datum=WGS84 '//&<br>   & '+ellps=WGS84 +towgs84=0,0,0'//char(0))<br><br>   SystemOutPtr = pj_init_plus('+proj=ob_tran +o_proj=eqc '//&<br>   & '+o_lat_p=30.0 +a=57.29578 +lon_0=0'//char(0))<br>!<br>! Read a point in lon-lat [degree], transform to [radians]<br>!<br>   READ(*,*) x,y<br>   WRITE(*,10) x,y<br> 10 FORMAT('In:  ',F15.5,1X,F15.5)<br><br>   z = 0._C_DOUBLE<br>   <br>   x = x * 0.0174533_C_DOUBLE<br>   y = y * 0.0174533_C_DOUBLE<br>!<br>! Transform to shifted-pole 60 coordinates<br>!<br>   MyResult = pj_transform(SystemInPtr,SystemOutPtr, 1, 0, x, y, z)<br>!<br>! Show results<br>!<br>   WRITE(*,20) x,y<br> 20 FORMAT('Out: ',F15.5,1X,F15.5)<br>!<br>   END PROGRAM TestProj<br><br><br>To compile it with g95:<br>g95 libmyp.f90 testproj.f90 -lproj -o testproj.exe<br><br><br>As sample input feed it:<br>
<font class="fixed_width" face="Courier, Monospaced">2.21    -6.44 <br> </font>
<p><font class="fixed_width" face="Courier, Monospaced">The output shold be:</font></p>
<p><font class="fixed_width" face="Courier, Monospaced">In:          3.70000        53.50000 <br> Out:         2.21382        -6.43806 <br> </font></p>
<br>
The program should now show the new coordinates after rotation of the South Pole<br>over 60 degrees northwards over the Greenwich meridian, which is reflected by the difference between the two y-coordinates.<br><br>Challenge for a bonus: the above solution uses ISO_C_BINDING. Not every compiler<br>supports this unit, which is an F2003 extension. Anyone who can show me how to do the same without ISO_C_BINDING is encouraged to do so!<br><br>Have fun!<br><br><br>Arjan<br><br></FONT>