１．はじめに

　この記事では、Ubuntuで Fortran から Proj4 の機能を使う方法をご紹介します。
　検証は、IntelのoneAPIで行いました。

２．Proj4 のインストール

２．１．Projとは

　Proj はGISソフトでも使用されている投影変換のライブラリです。proj-stringと呼ばれる文字列で投影方法を指定して、投影変換を行ってくれます。
Proj のサイト
　Proj は2025年4月現在、最新版はver9.6ですが、ver5からそれ以前とAPIが大きく変わりました。この記事では、使い方が簡単なver4以前のAPIを使用する方法をご紹介します。

２．２．Proj のパッケージ

　Ubuntuでは、libproj-dev という Proj の開発用パッケージが用意されていて、既にインストールされている環境が多いと思いますが、ver8以降ではver4のAPIが使えなくなっているため、ver4を別途インストールします。

２．３．oneAPI について

　事前にoneAPIをインストールし、icx（またはicc）と ifx（ifort）が実行できるようにしておいてください。Qiitaで「oneAPI ubuntu」で検索するといろいろ出てきます。

２．４．ソースのダウンロード

　ver4の最終版は2025年4月現在、ver4.9.3なので
Proj のダウンロードページ
から、proj-4.9.3.tar.gz をダウンロードしてきてください。

２．５．ビルドしてインストール

　ダウンロードした proj-4.9.3.tar.gz を展開可能な場所に置いて

tar xvzf proj-4.9.3.tar.gz
cd proj-4.9.3
./configure --prefix=$HOME/usr CC=icx
make
make install

とすると、$HOME/usrの下にインストールされます。icx は icc でも大丈夫です。
インストール先はどこでも良いですが、この記事では$HOME/usrにインストールした前提で説明します。

３．Fortran からの呼び出し

　インターフェース定義モジュール

proj4.f90

    module proj4

    use iso_c_binding
    implicit none

    private
    public      ::  pjf90_init_plus, pjf90_transform, pjf90_free, pjf90_is_latlong

    interface

    type(c_ptr) function pj_init_plus(cprojstring) bind(c, name='pj_init_plus')
        import c_char, c_ptr
        character(1, kind=c_char), intent(in) :: cprojstring(*)
    end function pj_init_plus

    integer(4) function pj_transform(prjsrc, prjdst, point_count, point_offset, x, y, z) bind(c, name='pj_transform')
        import c_ptr
        type(c_ptr), value, intent(in)  ::  prjsrc, prjdst, z
        integer(4), value, intent(in)   ::  point_count, point_offset
        real(8), intent(in)             ::  x, y
    end function pj_transform

    subroutine pjf90_free(prj) bind(c, name='pj_free')
        import c_ptr
        type(c_ptr), value, intent(in)  ::  prj
    end subroutine pjf90_free

    integer(4) function pjf90_is_latlong(prj) bind(c, name='pj_is_latlong')
        import c_ptr
        type(c_ptr), value, intent(in)  ::  prj
    end function pjf90_is_latlong

    end interface

    interface pjf90_transform
        module procedure pjf90_transformN, pjf90_transform1
    end interface pjf90_transform

    contains

    type(c_ptr) function pjf90_init_plus(cprojstring)
        implicit none
        character(len=*), intent(in)    ::  cprojstring
        pjf90_init_plus = pj_init_plus(trim(cprojstring)//c_null_char)
    end function pjf90_init_plus

    integer(4) function pjf90_transformN(prjsrc, prjdst, np, x, y)
        implicit none
        type(c_ptr), intent(in)         ::  prjsrc, prjdst
        integer(4), intent(in)          ::  np
        real(8), intent(inout)          ::  x(*), y(*)
        integer(4)          ::  ip
        real(8), parameter  ::  pi = 3.14159265358979323846d0
        if (pjf90_is_latlong(prjsrc) /= 0) then
            do ip=1, np
                x(ip) = x(ip) / 180d0 * pi
                y(ip) = y(ip) / 180d0 * pi
            end do
        end if
        pjf90_transformN = pj_transform(prjsrc, prjdst, np, 1, x(1), y(1), c_null_ptr)
        if (pjf90_is_latlong(prjdst) /= 0) then
            do ip=1, np
                x(ip) = x(ip) * 180d0 / pi
                y(ip) = y(ip) * 180d0 / pi
            end do
        end if
    end function pjf90_transformN

    integer(4) function pjf90_transform1(prjsrc, prjdst, x, y)
        implicit none
        type(c_ptr), intent(in)         ::  prjsrc, prjdst
        real(8), intent(inout)          ::  x, y
        real(8)     ::  xx(1), yy(1)
        xx(1) = x
        yy(1) = y
        pjf90_transform1 = pjf90_transformN(prjsrc, prjdst, 1, xx, yy)
        x = xx(1)
        y = yy(1)
    end function pjf90_transform1

    end module proj4

　オリジナルのpj_transformでは、Z座標が求められたり、地点データ間のオフセット量を指定できますが、それは使わないので除外しました。
　また、オリジナルのpj_transformでは、緯度・経度座標はラジアン値で渡す必要があるのですが、呼び出し元ではdegree値で管理したいので、pj_transformの呼び出しの前後でラジアン値とdegree値を変換すようにしました。

　テストプログラム

t_proj4.f90

    program t_proj4

    use iso_c_binding, only : c_ptr

    use proj4

    implicit none

    type(c_ptr)             ::  src, dst
    integer(4)              ::  ierr, ip
    integer(4), parameter   ::  np = 38
    real(8)                 ::  lons(np) = (/ 143.2586940d0, 143.2587500d0, 143.2587500d0, 141.6699440d0, 141.6699440d0, 141.0033330d0, 141.0033330d0, 140.4766940d0, 140.4766940d0, 140.2912780d0, 140.2912780d0, 141.1870830d0, 141.1870830d0, 139.9837220d0, 139.9837220d0, 139.8237220d0, 139.8237220d0, 140.5228890d0, 140.5228890d0, 140.3433610d0, 140.3433610d0, 141.4170280d0, 141.4170280d0, 141.9133060d0, 141.9133060d0, 143.7349440d0, 143.7349440d0, 144.7758060d0, 144.7758060d0, 145.3433330d0, 145.3433330d0, 145.2295280d0, 145.2295280d0, 145.8178610d0, 145.8178610d0, 143.9940830d0, 143.9940830d0, 143.2586940d0 /)
    real(8)                 ::  lats(np) = (/ 41.9991110d0, 41.9983060d0, 41.9983060d0, 42.6495280d0, 42.6495280d0, 42.3012220d0, 42.3012220d0, 42.5887220d0, 42.5887220d0, 42.2491670d0, 42.2491670d0, 41.8008060d0, 41.8008060d0, 41.5558060d0, 41.5558060d0, 42.6141670d0, 42.6141670d0, 42.9908610d0, 42.9908610d0, 43.3353890d0, 43.3353890d0, 43.3224720d0, 43.3224720d0, 45.5187220d0, 45.5187220d0, 44.1062500d0, 44.1062500d0, 43.9270830d0, 43.9270830d0, 44.3420560d0, 44.3420560d0, 43.3366670d0, 43.3366670d0, 43.3866390d0, 43.3866390d0, 42.9270560d0, 42.9270560d0, 41.9991110d0 /)
    real(8)                 ::  xx(np), yy(np)

    src = pjf90_init_plus('+proj=longlat +datum=WGS84')
    dst = pjf90_init_plus('+proj=tmerc +lat_0=44 +lon_0=142.25 +k=0.9999 +ellps=GRS80')

    xx = lons
    yy = lats

    ierr = pjf90_transform(src, dst, np, xx, yy)

    do ip=1, np
        write(*, '(2f20.2)') xx(ip), yy(ip)
    end do
    
    call pjf90_free(src)
    call pjf90_free(dst)

    stop

    contains

    end program t_proj4

　コンパイルはこの2つのファイルがあるフォルダで

ifx proj4.f90 t_proj4.f90 ${HOME}/usr/lib/libproj.a -o t_proj4.exe

とします。ifxの代わりにifortでも大丈夫です。
　この例では、北海道の輪郭のWGS84の緯度・経度座標値を、平面直角座標12系（JGD2000）値に変換します