[整理版]坐标转换从经纬度坐标到大地坐标及源码

[整理版]坐标转换从经纬度坐标到大地坐标及源码 坐标转换从经纬度坐标到大地坐标及源码 利用网络上开源的资料，可以很容易的实现从经纬度坐标向各种投影坐标的转 换，美国地质调查局开发USGS的GCTP就是很好的东西之一，有,语言版本的支 持各种投影类型的源代码，比如UNIVERSAL TRANSVERSE MERCATOR， ALBERS CONICAL EQUAL AREA ，LAMBERT CONFORMAL CONIC等等，我们机房使用的就是 LAMBERT CONFORMAL CONIC(LAMBERT 圆锥等角投影)，利用GCTP提供的源代码 详见source文件夹下的lamccfor.c 和 lamccinv.c ，可以轻松实现经纬度 投影坐标和LAMBERT CONFORMAL CONIC坐标的相互转化，其他投影方式使用的比 较少，还要进一步研究，就以后在写了( lamccforint函数设置LAMBERT CONFORMAL CONIC投影的各个参数，比如长半轴，短半轴，中心点经纬度坐标，标准纬线等等信息( long lamccforint(r_maj,r_min,lat1,lat2,c_lon,c_lat,false_east,false_north) double r_maj; /* majoraxis */ double r_min; /* minoraxis */ double lat1; /* first standard parallel */ double lat2; /* second standard parallel */ double c_lon; /* center longitude */ double c_lat; /* center latitude */ /* x offset in double false_east; meters */ double false_north; /* y offset in meters */ { double sin_po; /* sin value */ double cos_po; /* cos value */ double con; /* temporary variable */ double ms1; /* small m 1 */ double ms2; /* small m 2 */ double temp; /* temporary variable */ double ts0; /* small t 0 */ /* small t double ts1; 1 */ double ts2; /* small t 2 */ r_major = r_maj; r_minor = r_min; false_northing = false_north; false_easting = false_east; /* Standard Parallels cannot be equal and on opposite sides of the equator --------------------------------------------------------------------- ---*/ if (fabs(lat1 lat2) < EPSLN) { p_error("Equal latitudes for St. Parallels on opposite sides of equator", "lamcc-for"); return(41); } temp = r_minor / r_major; es = 1.0 - SQUARE(temp); e = sqrt(es); center_lon = c_lon; center_lat = c_lat; sincos(lat1,&sin_po,&cos_po); con = sin_po; ms1 = msfnz(e,sin_po,cos_po); ts1 = tsfnz(e,lat1,sin_po); sincos(lat2,&sin_po,&cos_po); ms2 = msfnz(e,sin_po,cos_po); ts2 = tsfnz(e,lat2,sin_po); sin_po = sin(center_lat); ts0 = tsfnz(e,center_lat,sin_po); if (fabs(lat1 - lat2) > EPSLN) ns = log (ms1/ms2)/ log (ts1/ts2); else ns = con; f0 = ms1 / (ns * pow(ts1,ns)); rh = r_major * f0 * pow(ts0,ns); /* Report parameters to the user -----------------------------*/ ptitle("LAMBERT CONFORMAL CONIC"); radius2(r_major, r_minor); stanparl(lat1,lat2); cenlonmer(center_lon); origin(c_lat); offsetp(false_easting,false_northing); return(OK); } lamccfor函数实现输入任意点经纬度值输出此投影下对应的点的大地坐标 long lamccfor(lon, lat, x, y) double lon; /* (I) Longitude */ double lat; /* (I) Latitude */ double *x; /* (O) X projection coordinate */ double *y; /* (O) Y projection coordinate */ { /* temporary angle double con; variable */ double rh1; /* height above ellipsoid */ double sinphi; /* sin value */ double theta; /* angle */ double ts; /* small value t */ con = fabs( fabs(lat) - HALF_PI); if (con > EPSLN) { sinphi = sin(lat); ts = tsfnz(e,lat,sinphi); rh1 = r_major * f0 * pow(ts,ns); } else { con = lat * ns; if (con <= 0) { p_error("Point can not be projected","lamcc-for"); return(44); } rh1 = 0; } theta = ns * adjust_lon(lon - center_lon); *x = rh1 * sin(theta) false_easting; *y = rh - rh1 * cos(theta) false_northing; return(OK); } lamccinv函数的作用与lamccfor函数相反，计算从大地坐标变换到对应的经纬度值 long lamccinv(x , y, lon, lat) double x; /* (O) X projection coordinate */ double y; /* (O) Y projection coordinate */ double *lon; /* (I) Longitude */ double *lat; /* (I) Latitude */ { double rh1; /* height above ellipsoid */ double con; /* sign variable */ double ts; /* small t */ double theta; /* angle */ long flag; /* error flag */ flag = 0; x -= false_easting; y = rh - y false_northing; if (ns > 0) { rh1 = sqrt (x * x y * y); con = 1.0; } else { rh1 = -sqrt (x * x y * y); con = -1.0; } theta = 0.0; if (rh1 != 0) theta = atan2((con * x),(con * y)); if ((rh1 != 0) || (ns > 0.0)) { con = 1.0/ns; ts = pow((rh1/(r_major * f0)),con); *lat = phi2z(e,ts,&flag); if (flag != 0) return(flag); } else *lat = -HALF_PI; *lon = adjust_lon(theta/ns center_lon); return(OK); }