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
如下圖所示，本文已經包含了橢球的七參數和三參數變換（不同基準的變換）和投影變換。
