緯度経度⇔投影座標値との相互変換を行うライブラリです。
ほぼ世界中の座標系と、マニアックな投影法までをカバーしてあります。
GRASSやMapServer等のオープンソースのGISソフトでは必ずといっていいほど使用されています。

コンパイル

ver. 4.5.0

./fonfigure → makeで普通にできますが、DLLが作れません。
nmakeで作ることはできるようです。この場合はsrcフォルダに移動してnmake -f makefile.vcを実行してください。

準備

必要なヘッダファイルはprojects.hとproj_api.hで、projects.hをインクルードします。
同時にwindows.hをインクルードしていると、PVALUE型が衝突するので、適当な名前に変更しておきます。

pj_init()関数

まず、pj_init()関数によって投影法、楕円体、その投影座標での原点の緯度経度を指定します。
指定方法の例は、以下のとおり。

static char *params[]={
  "proj=utm",
  "ellps=GRS80",
  "lon_0=135",
  "no_defs"
};

PJ *ref;

if ( !( ref = pj_init( sizeof(params)/sizeof((char *), params) ) )
{
  fprintf( stderr, "%sn", pj_strerrno( pj_errno ) );
}

またはpj_init_plus()関数を使って以下のように設定することもできます。

  sprintf( szInit,
  "+proj=tmerc +lat_0=%s +lon_0=%s +k=0.999900 +x_0=0 +y_0=0 
    +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
    szLat[iSystem], szLon[iSystem] );

  if ( !( ref = pj_init_plus( szInit ) ) )
  {
    ::MessageBox( NULL, "座標系の初期化に失敗しました", "Error", MB_OK );
    return;
  }

1行目は投影方法、2行目は楕円体、3,4行目は投影座標の中心を指定するらしい。
特に指定しなくてもいい場合があるらしく、その場合はno_defsを指定します。
たとえば、UTM座標値を求めたい場合などは、北半球では緯度の設定は不要です。

各項目で、設定することができる値は以下を参照。

投影法一覧
楕円体一覧

proj4のソースコードの中にはnadフォルダの中にepsgというファイルがあります。
この中身は各座標系のepsgコードとそのパラメータが記述されています。
なので、MapServerのmapファイルの記述のようにpj_init_plus( "epsg=4724" );などと指定できそうなものですが、どうもできないようです。
EPSGコードで指定したい場合はpj_init_plus("+init=epsg:4612");とします。

座標値の求め方

pj_init()関数が成功したら、ライブラリはその投影法と楕円体、原点を元に変換するようになります。
緯度経度から座標値を求めるには、以下のようにする。

projUV data;

data.v *= DEG_TO_RAD;
data.u *= DEG_TO_RAD;
data = pj_fwd( data, ref );
if ( data.u != HUGE_VAL )
{
	printf( "%.3ft%.3fn", data.u, data.v );
}
else
{
	fprintf( stderr, "!2n" );
}

gdalを使うなら

gdalを使うと、EPSG番号で変換することができます。
内部的にはgdalからproj4を呼んでいます。

この場合のコードは大体以下のようにします。
gdalがエクスポートしているクラスを使ってもいいはずですが、
筆者の環境ではなぜかOGRSpatialReferenceオブジェクトをdelete
するときに謎の赤バツが出てしまいます。
（何かしでかしているんだとは思いますが、、、)

#include "ogr_srs_api.h"

int main( int argc, char *argv[] )
{
	OGRSpatialReferenceH hIn, hOut;
	OGRCoordinateTransformationH hCT;
	
	double aX[] = { 0.0, 1000.0, -1000.0, 2000.0 };
	double aY[] = { 0.0, 1000.0, -1000.0, 2000.0 };

	hIn = OSRNewSpatialReference( NULL );
	hOut = OSRNewSpatialReference( NULL );

	OSRImportFromEPSG( hIn, 2451 );
	OSRImportFromEPSG( hOut, 4319 );

	hCT = OCTNewCoordinateTransformation( hIn, hOut );

	OCTTransform( hCT, 4, aX, aY, NULL );
	
	for ( int i = 0; i < 4; i++ )
	{
		printf( "%8.3f, %8.3fn", aX[i], aY[i] );
	}

	OCTDestroyCoordinateTransformation( hCT );

	OSRDestroySpatialReference( hIn );
	OSRDestroySpatialReference( hOut );
	
	return 0;
}