<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>