GIS坐标转换库Proj.4的使用

     GIS(Geographic Information System 地理信息系统）领域中最常提及 的一个概念是坐标系统，当我们提及一个地理位置的时候，与之伴随而产生的是该位置必定在一个空间参考下。当我们使用GPS设备获取到某个位置的经纬度的时候，我们实际上得到的是一个WGS-84椭球坐标系统下的空间地理坐标。目前空间坐标系统已经十分繁多，这些坐标系统有着特定的使用场合，并且可以通过一定的方式实现它们之间的相互转换，本文提到的Proj.4库正是这样一个提供了不同坐标系统（空间参考）之间相互转换的强大的工具。关于空间参考的详细介绍以及它们之间的转换方法可以参考以下的内容：Geometric Aspects of Mapping

• Proj.4库的安装

本文主要是介绍Proj.4库的使用，首先先安装Proj.4，可以去Porj.4的官网下载源码，目前Proj.4的源码已经迁移到Github，具体地址是： https://github.com/OSGeo/proj.4，代码下载之后得到的目录结构如下图所示：

下载之后就可以开始编译源码了，可以使用目录中的CMakeLists.txt，使用cmake的工具进行生成并编译。也可以按照目录中的README文件中提到的方式在命令行中进行编译，本文采用后者，因为这种方式非常方便就可以生成可供开发的SDK。如果读者需要去了解Proj.4的实现并需要亲自去修改Proj.4的源码，那么使用前一种方式将Proj.4组织成Visual Studio或者XCode等IDE的解决方案和工程是十分必要的。编译生成方式如下：（本文在Windows下进行，使用Visual Studio 2013来进行编译）

1. 首先在开始菜单中找到Visual Studio 2013的命令行工具，并启动它：

2. 运行下面的命令

C:\> cd proj   //进入Proj.4源码目录（以自己的位置为准）
C:\PROJ> nmake /f makefile.vc
C:\PROJ> nmake /f makefile.vc install-all

3.由于使用的是默认的命令，因此最终编译的结果会出现在C盘根目录下的PROJ目录下，目录结构如下：

至此编译过程完成。可以查看各个目录中的内容： bin目录是proj库生成的一些使用工具以及动态链接库，include是我们进行开发需要的头文件，lib则是链接库（包含静态链接和动态链接中的符号库），share目录中提供了一些数据，包含常用的坐标系统定义等

• Proj.4工具的使用

在Proj4编译之后，在bin目录中会生成一些工具，包括如下几个工具：

1. cs2cs.exe

2. geod.exe

3. nad2bin.exe

4. proj.exe

下面我们来看一下这些工具都有什么作用以及如何使用它们：

• cs2cs.exe

cs2cs工具可以实现两个坐标系统之间的互相转换，这两个坐标系统可以是投影坐标与地理坐标之间的转换，也可以是基准面之间的转换，cs2cs的参数格式如下：

cs2cs [-eEfIlrstvwW [args] ] [+opts[=args] ] [ +to [+opts [=arg] ] ] file[s]

• proj.exe

proj工具用来转换地理坐标到投影坐标以及投影坐标到地理坐标，但是有一个前提条件是这两个坐标系统在同一个基准面之下。 参考： http://geoffair.net/projects/proj4.htm

• geod.exe

geod工具可以计算大圆的距离

• nad2bin

nad2bin可以将坐标投影信息转换为二进制的文件

• 使用Proj.4库

使用Proj.4库就像使用常用的C或者C++库一样，在Visual Studio下创建工程，并把头文件和链接文件引入就可以进行开发了。使用Proj4进行开发，我们对以下的一些转换提供支持：

• 相同基准面之间的转换：投影坐标到地理坐标（经纬度）的转换以及反向转换（从地理坐标到投影坐标之间的转换）

这个功能实际上和Proj.4提供的proj.exe的功能是一样的，我们看看如何做到这一点。Proj4提供的开发接口并不多，但是接口参数中的坐标系统却十分繁多，假设我们需要实现将WGS84经纬度的坐标转换到墨卡托投影坐标系统下，代码如下：

#include <proj_api.h>
#include <iostream>


int main(int argc, char **argv) 
{
	projPJ pj_merc, pj_latlong;
	double x, y;

	if (!(pj_merc = pj_init_plus("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")))
		exit(1);
	if (!(pj_latlong = pj_init_plus("+proj=longlat +datum=WGS84 +no_defs")))
		exit(1);

	x = -9.866554;
	y = 7.454779;

	x *= DEG_TO_RAD;
	y *= DEG_TO_RAD;

	pj_transform(pj_latlong, pj_merc, 1, 1, &x, &y, NULL);

	std::cout.precision(12);
	std::cout << "(" << x << " , " << y << ")" << std::endl;
	exit(0);
}

基本上我们需要调用的是proj4提供的函数 pj_transform函数，它的函数原型如下：

int pj_transform( projPJ src, projPJ dst, long point_count, int point_offset,
                  double *x, double *y, double *z );

我们传入 需要转换数据的坐标系统（源）和转换到指定的坐标系统（目标），由于上面的示例程序中只有一个点需要转换，并且这些数据之间也不存在着间隔，都是紧密排列的，因此后面的两个参数设置为1， x, y, z指针作为传出参数，将最后转换的结果写到x，y，z之中。

当我们进行反向转换的时候，只需要将源和目标的参数对换就可以了（需要注意有部分变换是不能反向转换的)

#include <proj_api.h>
#include <iostream>

//(-1098339.76716 , 826673.618947)
//	x = -9.866554;
//y = 7.454779;

int main(int argc, char **argv) 
{
	projPJ pj_merc, pj_latlong;
	double x, y;

	if (!(pj_merc = pj_init_plus("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")))
		exit(1);
	if (!(pj_latlong = pj_init_plus("+proj=longlat +datum=WGS84 +no_defs")))
		exit(1);

	x = -1098339.76716;
	y = 826673.618947;

	pj_transform(pj_merc, pj_latlong, 1, 1, &x, &y, NULL);

	std::cout.precision(12);
	std::cout << "(" << x /DEG_TO_RAD<< " , " << y/DEG_TO_RAD << ")" << std::endl;
	exit(0);
}

• 基准面之间的变换

在实现基准面之间的变换的时候，可以有两种方式：

1. 使用地心坐标系统进行变换，这就要求首先必须将地理坐标转换为地心坐标，由于地心坐标是三维的直角坐标系统，这样两个直角坐标之间的变换就可以使用3参数或者7参数进行变换

2. 直接使用地理坐标进行基准面变换，这种变换方式有相应的算法实现，目前这种方式据我个人目前所知，在Proj4中并没有给出实现

下面看看在Proj4里面怎么实现基准面之间的变换（第一种实现方式），可以使用Proj4提供的函数

int pj_datum_transform( projPJ src, projPJ dst, long point_count, int point_offset,
                        double *x, double *y, double *z );

x, y, z的输入和输出均是 BLH(经纬高程)，src和dst的参数必须包含基准面转换的参数（+towgs84)：可以是三参数或者七参数，类似于：

+towgs84=-199.87,74.79,246.62
+towgs=0,0,4.5,0,0,0.554,0.219

• 球面坐标到地心坐标的转换以及反转换

当我们进行基准面的三参数或者七参数计算的时候，需要将球面的地理坐标转换为地心三维直角坐标系下的坐标，在Proj4中给出了下面的函数：

int pj_geocentric_to_geodetic( double a, double es,
                               long point_count, int point_offset,
                               double *x, double *y, double *z );
int pj_geodetic_to_geocentric( double a, double es,
                               long point_count, int point_offset,
                               double *x, double *y, double *z );

在转换为地心坐标之后，我们可以根据两个不同地心坐标系下面得到几个已知控制点计算出它们之间的三参数（需要1个已知控制点对）或者七参数（需要三个已知控制点对）

基本上Proj4主要提供的转换函数就是这些，除此之外它提供了初始化、析构以及判断坐标系是投影或者球面地理坐标系的函数。下面我们给出两个不同基准面的投影坐标之间的转换（如北京54坐标到西安80坐标之间的转换）过程，如下图所示：