proj4投影库的使用与安装

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

• 使用Proj.4库

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

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

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

[cpp] view
plain copy

1. #include <proj_api.h>  
2. #include <iostream>  
3.   
4.   
5. int main(int argc, char **argv)   
6. {  
7.     projPJ pj_merc, pj_latlong;  
8.     double x, y;  
9.   
10.     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”)))  
11.         exit(1);  
12.     if (!(pj_latlong = pj_init_plus(“+proj=longlat +datum=WGS84 +no_defs”)))  
13.         exit(1);  
14.   
15.     x = -9.866554;  
16.     y = 7.454779;  
17.   
18.     x *= DEG_TO_RAD;  
19.     y *= DEG_TO_RAD;  
20.   
21.     pj_transform(pj_latlong, pj_merc, 1, 1, &x, &y, NULL);  
22.   
23.     std::cout.precision(12);  
24.     std::cout << “(“ << x << ” , “ << y << “)” << std::endl;  
25.     exit(0);  
26. }  

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

[cpp] view
plain copy

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

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

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

[cpp] view
plain copy

1. #include <proj_api.h>  
2. #include <iostream>  
3.   
4. //(-1098339.76716 , 826673.618947)  
5. //  x = -9.866554;  
6. //y = 7.454779;  
7.   
8. int main(int argc, char **argv)   
9. {  
10.     projPJ pj_merc, pj_latlong;  
11.     double x, y;  
12.   
13.     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”)))  
14.         exit(1);  
15.     if (!(pj_latlong = pj_init_plus(“+proj=longlat +datum=WGS84 +no_defs”)))  
16.         exit(1);  
17.   
18.     x = -1098339.76716;  
19.     y = 826673.618947;  
20.   
21.     pj_transform(pj_merc, pj_latlong, 1, 1, &x, &y, NULL);  
22.   
23.     std::cout.precision(12);  
24.     std::cout << “(“ << x /DEG_TO_RAD<< ” , “ << y/DEG_TO_RAD << “)” << std::endl;  
25.     exit(0);  
26. }  

• 基准面之间的变换

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

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

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

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

[cpp] view
plain copy

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

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

[cpp] view
plain copy

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

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

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

[cpp] view
plain copy

1. int pj_geocentric_to_geodetic( double a, double es,  
2.                                long point_count, int point_offset,  
3.                                double *x, double *y, double *z );  
4. int pj_geodetic_to_geocentric( double a, double es,  
5.                                long point_count, int point_offset,  
6.                                double *x, double *y, double *z );  

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

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

转载自：https://blog.csdn.net/jiangbingbo123/article/details/51062349