１．はじめに
この記事では、Ubuntuで Intel Fortran から Proj4 の機能を使う方法をご紹介します。
ほんとは iso_c_binding など使って gfortran でも使えるサンプルを用意したかったんですが、勉強不足なため使い慣れた Intel Fortran のディレクティブで実装しました。
２．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が実行できるようにしておいてください。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にインストールした前提で説明します。
３．Intel Fortran からの呼び出し
インターフェース定義モジュール
module proj4
implicit none
private
public :: pjf90_init_plus, pjf90_transform, pjf90_free, pjf90_is_latlong
interface
integer(8) function pj_init_plus(cprojstring)
!dec$ attributes c,reference,alias:'pj_init_plus' :: pj_init_plus
character(len=*), intent(in) :: cprojstring
!dec$ attributes reference :: cprojstring
end function pj_init_plus
integer(4) function pj_transform(src, dst, point_count, point_offset, x, y, z)
!dec$ attributes c,reference,alias:'pj_transform' :: pj_transform
integer(8), intent(in) :: src, dst, z
integer(4), intent(in) :: point_count, point_offset
real(8), intent(in) :: x, y
!dec$ attributes value :: src, dst, point_count, point_offset, z
!dec$ attributes reference :: x, y
end function pj_transform
subroutine pjf90_free(prj)
!dec$ attributes c,reference,alias:'pj_free' :: pjf90_free
integer(8), intent(in) :: prj
!dec$ attributes value :: prj
end subroutine pjf90_free
integer(4) function pjf90_is_latlong(prj)
!dec$ attributes c,reference,alias:'pj_is_latlong' :: pjf90_is_latlong
integer(8), intent(in) :: prj
!dec$ attributes value :: prj
end function pjf90_is_latlong
end interface
interface pjf90_transform
module procedure pjf90_transformN, pjf90_transform1
end interface pjf90_transform
contains
integer(8) function pjf90_init_plus(cprojstring)
implicit none
character(len=*), intent(in) :: cprojstring
pjf90_init_plus = pj_init_plus(trim(cprojstring)//char(0))
end function pjf90_init_plus
integer(4) function pjf90_transformN(src, dst, np, x, y)
implicit none
integer(8), intent(in) :: src, dst
integer(4), intent(in) :: np
real(8), intent(inout) :: x(*), y(*)
integer(4) :: ip
real(8), parameter :: pi = 3.14159265358979323846d0
if (pjf90_is_latlong(src) /= 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(src, dst, np, 1, x(1), y(1), 0)
if (pjf90_is_latlong(dst) /= 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(src, dst, x, y)
implicit none
integer(8), intent(in) :: src, dst
real(8), intent(inout) :: x, y
real(8) :: xx(1), yy(1)
xx(1) = x
yy(1) = y
pjf90_transform1 = pjf90_transformN(src, dst, 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値を変換すようにしました。
テストプログラム
program t_proj4
use proj4
implicit none
integer(8) :: 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）値に変換します