坐标转换从经纬度坐标到大地坐标及源码.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)

   {

   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  */

{

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)

    {

    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);

 }