proj4测绘坐标系转换
相同基准面之间投影坐标和地理坐标的转换
pj_transform
projPJ pj_merc,pj_latlong;
double x,y;
//定义墨卡托投影坐标系
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");
//定义地理坐标系
pj_latlong = pj_init_plus("+proj=longlat +datum=WGS84 +no_defs");
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);
cout.precision(12);
cout << "(" << x << " , " << y << ")" << endl;
//投影坐标转地理坐标
pj_transform(pj_merc,pj_latlong,1,1,&x,&y,NULL);
x*=RAD_TO_DEG;
y*=RAD_TO_DEG;
cout.precision(8);
cout << "(" << x << " , " << y << ")" << endl;
相同基准下地理坐标和空间直角坐标系的转换（即大地坐标与地心坐标的转换）
pj_geocentric_to_geodetic
：地心坐标转大地坐标（XYZ->BLH）
pj_geodetic_to_geocentric
：大地坐标转地心坐标（BLH->XYZ）
WGS84基准下的转换
typedef struct { double u, v,w; } projUVW;
double wgs84_a = 6378137.0;//长半轴
double wgs84_f = 1.0/298.257223563;//扁率
double wgs84_b = wgs84_a * (1 - wgs84_f);//短半轴
double wgs84_e=(wgs84_a*wgs84_a-wgs84_b*wgs84_b)/(wgs84_a*wgs84_a);//第一偏心率的平方
projUVW pt3[2] = {{116.5164884833, 39.7663036028, 23.220}, {116.987654321, 39.453, 18}};
for (int i = 0; i < 2; i++)
{
cout << "第" << i << "个点的转换：\n";
pt3[i].u *= DEG_TO_RAD;
pt3[i].v *= DEG_TO_RAD;
pj_geodetic_to_geocentric(wgs84_a, wgs84_e, 1, 1, &pt3[i].u, &pt3[i].v, &pt3[i].w);
cout.setf(ios::fixed);
cout.precision(4);
cout.width(12);
cout << "\t经纬度转换为空间坐标\n\t";
cout << pt3[i].u << "\t" << pt3[i].v << "\t" << pt3[i].w << endl;
pj_geocentric_to_geodetic(wgs84_a, wgs84_e, 1, 1, &pt3[i].u, &pt3[i].v, &pt3[i].w);
cout.setf(ios::fixed);
cout.precision(9);
cout.width(12);
cout << "\t空间坐标转换为经纬度\n\t";
cout << pt3[i].u * RAD_TO_DEG << "\t" << pt3[i].v * RAD_TO_DEG << "\t" ;
cout.precision(4);
cout << pt3[i].w << endl;
}
等价于：
typedef struct { double u, v,w; } projUVW;
// 地心坐标系
const char* geoccs="+proj=geocent +datum=WGS84";
// 大地坐标系
const char* latlon="+proj=latlong +datum=WGS84";
projPJ pjGeoccs, pjLatlon;
//初始化当前投影对象
pjGeoccs= pj_init_plus(geoccs);
pjLatlon= pj_init_plus(latlon);
projUVW pt3[2] = {{116.5164884833, 39.7663036028, 23.220}, {116.987654321, 39.453, 18}};
for (int i = 0; i < 2; i++)
{
cout << "第" << i << "个点的转换：\n";
pt3[i].u *= DEG_TO_RAD;
pt3[i].v *= DEG_TO_RAD;
pj_transform(pjLatlon,pjGeoccs, 1, 1, &pt3[i].u, &pt3[i].v, &pt3[i].w);
cout.setf(ios::fixed);
cout.precision(4);
cout.width(12);
cout << "\t经纬度转换为空间坐标\n\t";
cout << pt3[i].u << "\t" << pt3[i].v << "\t" << pt3[i].w << endl;
pj_transform(pjGeoccs, pjLatlon, 1, 1, &pt3[i].u, &pt3[i].v, &pt3[i].w);
cout.setf(ios::fixed);
cout.precision(9);
cout.width(12);
cout << "\t空间坐标转换为经纬度\n\t";
cout << pt3[i].u * RAD_TO_DEG << "\t" << pt3[i].v * RAD_TO_DEG << "\t" ;
cout.precision(4);
cout << pt3[i].w << endl;
}
不同基准下空间直角坐标系下的转换
pj_datum_transform
三参数
- GRS80下的地理坐标转GRS80下的空间直角坐标系
- 三参数将GRS80空间直角坐标系转WGS84空间直角坐标系
- WGS84空间直角坐标系转WGS84地理坐标系
//grs椭球参数
double grs80_a=6378137.0;WGS
double grs80_f=1.0/298.257222101;
double grs80_b=grs80_a*(1-grs80_f);
double grs80_e=(grs80_a*grs80_a-grs80_b*grs80_b)/(grs80_a*grs80_a) ;
//BLH
x=20;
y=35;
z=1000;
x*=DEG_TO_RAD;
y*=DEG_TO_RAD;
pj_geodetic_to_geocentric(grs80_a, grs80_e, 1, 1, &x, &y, &z);
//三参数 xyz轴的偏移
x-=199.87;
y+=74.79;
z+=246.62;
//wgs84椭球参数
double wgs84_a = 6378137.0;
double wgs84_f = 1.0/298.257223563;
double wgs84_b = wgs84_a * (1 - wgs84_f);
double wgs84_e=(wgs84_a*wgs84_a-wgs84_b*wgs84_b)/(wgs84_a*wgs84_a);
pj_geocentric_to_geodetic(wgs84_a,wgs84_e,1,1,&x,&y,&z);
x*=RAD_TO_DEG;
y*=RAD_TO_DEG;
cout.precision(12);
cout << "(" << x << " , " << y << ")" << endl;
等价于：
+towgs84=-199.87,74.79,246.62
projPJ pj_grs80,pj_wgs84;
pj_grs80=pj_init_plus("+proj=latlong +ellps=GRS80 +towgs84=-199.87,74.79,246.62");
pj_wgs84 = pj_init_plus("+proj=latlong +datum=WGS84");
double x,y,z;
x=20;
y=35;
z=1000;
x*=DEG_TO_RAD;
y*=DEG_TO_RAD;
pj_transform(pj_grs80,pj_wgs84,1,1,&x,&y,&z);
x*=RAD_TO_DEG;
y*=RAD_TO_DEG;
cout.precision(12);
cout << "(" << x << " , " << y << ")" << endl;
七参数
- WGS72下的地理坐标转WGS72下的空间直角坐标系
- 七参数将WGS72空间直角坐标系转WGS84空间直角坐标系
- WGS84空间直角坐标系转WGS84地理坐标系
#define PI 3.14159265358979323e0
#define SEC_TO_RAD PI/180/3600
//WGS72 椭球参数
double wgs72_f=1.0/298.26;
double wgs72_b=wgs72_a*(1-wgs72_f);
double wgs72_e=(wgs72_a*wgs72_a-wgs72_b*wgs72_b)/(wgs72_a*wgs72_a) ;//第一偏心率的平方
x=4;
y=55;
z=1000;
x*=DEG_TO_RAD;
y*=DEG_TO_RAD;
pj_geodetic_to_geocentric(wgs72_a, wgs72_e, 1, 1, &x, &y, &z);
// 七参数
double Dx_BF =0;
double Dy_BF =0;
double Dz_BF =4.5;
double Rx_BF =0;
double Ry_BF =0;
double Rz_BF =0.554*SEC_TO_RAD;//ms
double M_BF =1.0+0.219/1000000.0;//尺度变化值(实际缩放值=1.0+缩放因子/1000000)
double x_out = M_BF*( x - Rz_BF*y + Ry_BF*z) + Dx_BF;
double y_out = M_BF*( Rz_BF*x + y - Rx_BF*z) + Dy_BF;
double z_out = M_BF*(-Ry_BF*x + Rx_BF*y + z) + Dz_BF;
//wgs84椭球参数
double wgs84_a = 6378137.0;d
double wgs84_f = 1.0/298.257223563;
double wgs84_b = wgs84_a * (1 - wgs84_f);
double wgs84_e=(wgs84_a*wgs84_a-wgs84_b*wgs84_b)/(wgs84_a*wgs84_a);
pj_geocentric_to_geodetic(wgs84_a,wgs84_e,1,1,&x_out,&y_out,&z_out);
x_out*=RAD_TO_DEG;
y_out*=RAD_TO_DEG;
cout.precision(12);
cout << "(" << x_out << " , " << y_out << ")" << endl;
等价于：
+towgs84=0,0,4.5,0,0,0.554,0.219
projPJ pj_wgs72,pj_wgs84;
pj_wgs72=pj_init_plus("+proj=latlong +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.219");
pj_wgs84 = pj_init_plus("+proj=latlong +datum=WGS84");
x=4;
y=55;
z=1000;
x*=DEG_TO_RAD;
y*=DEG_TO_RAD;
pj_transform(pj_wgs72,pj_wgs84,1,1,&x,&y,&z);
亦等价于：
pj_datum_transform
projPJ pj_wgs72,pj_wgs84;
pj_wgs72=pj_init_plus("+proj=latlong +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.219");
pj_wgs84 = pj_init_plus("+proj=latlong +datum=WGS84");
x=4;
y=55;
z=1000;
x*=DEG_TO_RAD;
y*=DEG_TO_RAD;
pj_datum_transform(pj_wgs72,pj_wgs84,1,1,&x,&y,&z);
总结
坐标系转换的流程大抵如此：
pj_transform
->
pj_geodetic_to_geocentric
->
pj_datum_transform
->
pj_geocentric_to_geodetic
->
pj_transform
如下图所示，本文已经包含了椭球的七参数和三参数变换（不同基准的变换）和投影变换。
