坐标转换从经纬度坐标到大地坐标及源码.doc

坐标转换从经纬度坐标到大地坐标及源码利用网络上开源的资料，可以很容易的实现从经纬度坐标向各种投影坐标的转换，美国地质调查局开发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 */double false_east; /* x offset in 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 */double ts1; /* small t 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) 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 */double con; /* temporary angle 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) 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