Proj.NETProj.NET简介和示例(转载)

 Proj.NET ( http://www.codeplex.com/ProjNET )是一个.NET下开源的空间参照和投影引擎，遵循OGC相关标准。负责人(Coordinators )是 D_Guidi 和 SharpGIS，开发者(Developers)还有codekaizen、rstuven等，这些人也是.NET下其他的开源系统如GeoAPI.NET、SharpMap等主要贡献者。

 

    Proj.NET支持基准面转换，地理坐标系，投影坐标系，地心坐标系，可在多种.NET框架下使用，并支持SilverLight。它可进行点对点之间的坐标变换，同时也可以把坐标系转换为 Well-Known Text (WKT) 和 XML。目前支持投影类型有：Mercator、横轴墨卡托投影(Transverse Mercator)、Albers、正轴等角割圆锥投影(Lambert Conformal Conic、兰勃特投影）等。

 

    目前最新版发布版本是1.2，2009年5月，binaries内含Proj.NET.dll和Proj.NET.XML两个文件，创建程序时引入dll文件，即可使用相关类、函数、接口等。若想深入了解，建议把源码和示例一并下载。下面以实例来注解Proj.NET进行坐标转换的用法。在使用 Proj.NET 之前最后对地图投影基本概念，如投影参数，坐标系，SRID，WKT等概念有一定的了解。

 

下面有4对坐标要进行转换：

投影参数：椭球体Krasovsky_1940；Datum：北京1954；投影：兰勃特双标准纬线，25N，47N。
坐标数据： 20634500,4660000； 20635500,466000； 20634500, 4659000；20635500，4659000 。

 注：

    求经纬度无非就是从自定义投影坐标系转换到WGS1984坐标系，

    反之若求地理坐标则是从WGS84geogCS转换到要求的投影坐标系。

 

 首先了解下在Proj.NET创建坐标系的方法：

（1）Predefined （2）SRID （3）WKT （4）Code

·Predefined

  用预定义方法或SRID创建WGS84地理坐标系（预定义方法，十分简单，但系统预定义坐标系为数不多）：
  GeographicCoordinateSystem toCS = GeographicCoordinateSystem.WGS84;
·SRID
  使用SRID也比较容易，前提是引入EPSG的数据集，还有SridReader类文件，在源码中有

  ICoordinateSystem toCS = SRIDReader.GetCSbyID(4326);

·WKT

  用WKT创建给定的投影坐标系(写代码时注意转义字符"和删除换行)：
  stringbj1954Lcc = PROJCS["liongg",
      GEOGCS["GCS_Beijing_1954",
          DATUM["D_Beijing_1954",SPHEROID["Krasovsky_1940",6378245.0,298.3]],
          PRIMEM["Greenwich",0.0],
          UNIT["Degree",0.0174532925199433]],
      PROJECTION["Lambert_Conformal_Conic"],
      PARAMETER["False_Easting",20500000.0],
      PARAMETER["False_Northing",0.0],
      PARAMETER["Central_Meridian",0.0],
      PARAMETER["Standard_Parallel_1",25.0],
      PARAMETER["Standard_Parallel_2",47.0],
      PARAMETER["Scale_Factor",1.0],
      PARAMETER["Latitude_Of_Origin",0.0],
      UNIT["Meter",1.0]] ;
  IProjectedCoordinateSystem fromCS = CoordinateSystemWktReader.Parse(bj1954Lcc) as IProjectedCoordinateSystem;

·Code

  所谓Code方式即把WKT中的“元素”如DATUM、SPHEROID、UNIT等分别定义，最后组合。

  CoordinateSystemFactory cFac = new CoordinateSystemFactory();

  IEllipsoid ellipsoid = cFac.CreateFlattenedSphere("Bessel 1840", 6377397.155, 299.15281, LinearUnit.Metre);

  IHorizontalDatum datum = cFac.CreateHorizontalDatum("Bessel 1840", DatumType.HD_Geocentric, ellipsoid, null);

  IGeographicCoordinateSystem gcs = cFac.CreateGeographicCoordinateSystem("Bessel 1840", AngularUnit.Degrees, 

      datum, PrimeMeridian.Greenwich, new AxisInfo("Lon", AxisOrientationEnum.East),

      new AxisInfo("Lat", AxisOrientationEnum.North));

  List<ProjectionParameter> parameters = new List<ProjectionParameter>(5);

  parameters.Add(new ProjectionParameter("latitude_of_origin", 0));

  parameters.Add(new ProjectionParameter("central_meridian", 110));

  parameters.Add(new ProjectionParameter("scale_factor", 0.997));

  parameters.Add(new ProjectionParameter("false_easting", 3900000));

  parameters.Add(new ProjectionParameter("false_northing", 900000));

  IProjection projection = cFac.CreateProjection("Mercator_1SP", "Mercator_1SP", parameters);

  IProjectedCoordinateSystem coordsys = cFac.CreateProjectedCoordinateSystem("Makassar / NEIEZ", gcs,

      projection, LinearUnit.Metre, new AxisInfo("East", AxisOrientationEnum.East),

      new AxisInfo("North", AxisOrientationEnum.North));

 

坐标转换函数

转换多个坐标对的函数原型：List<double[]> TransformList(List<double[]> points);
参数和返回值都是List型，List的每个元素是一个一维double型数组。每个一维数组即一个坐标对：两个或三个元素，取决于有没有高程等第三个坐标Z，顺序为(x,y,z)。因此我们需要把给的坐标存储到List中。

 

下面是一个多点到多点的转换函数：

参数分别是源坐标系，目标坐标系，要转换的坐标值List，返回的坐标值List。

static void PtsToPts(ICoordinateSystem fromCS, ICoordinateSystem toCS,

    List<double[]> pts, out List<double[]> results)

{

    CoordinateTransformationFactory ctfac = new CoordinateTransformationFactory();

    ICoordinateTransformation trans = ctfac.CreateFromCoordinateSystems(fromCS, toCS);

    results = trans.MathTransform.TransformList(pts);

}

 

下面是主函数，简略起见，命名空间的引用，异常捕捉等省略，并不再注释。

static void Main(string[] args)

{

    string bj1954Lcc= "……"; //见上文

    IProjectedCoordinateSystem fromCS = CoordinateSystemWktReader.Parse(bj1954Lcc) as IProjectedCoordinateSystem;

    GeographicCoordinateSystem toCS = GeographicCoordinateSystem.WGS84;

    double[] xPts = new double[4] { 20634500, 20635500, 20634500, 20635500 };

    double[] yPts = new double[4] { 4660000, 4660000, 4659000, 4659000 };

    List<double[]> pts = new List<double[]>();

    for (int i = 0; i < 4; i++)

    {

        double[] xy = new double[] { xPts[i], yPts[i] };

        pts.Add(xy);

        List<double[]> results = new List<double[]>();

        PtsToPts(fromCS, toCS, pts, out results);

        foreach (double[] result in results)

        {

            Console.WriteLine("{0},{1}", result[0], result[1]);

        }

    }

}

下面是计算结果，可以看出来和使用proj.4的计算结果只有精度位数的差异。
118.611065451529,40.4621006468135
118.623041399559,40.4619479215536
118.610866262964,40.4529532936788
118.622840730867,40.4528005852267

 

proj.4是用C写的，很好很强大，我暂且只会简单使用，没看过源代码。Proj.4的计算结果：
Beijing 1954:   20634500.0   4660000.00
pj_inv result: 118.611065      40.462077

Beijing 1954:   20635500.0   4660000.00
pj_inv result: 118.623041      40.461924

Beijing 1954:   20634500.0   4659000.00
pj_inv result: 118.610866      40.452929

Beijing 1954:   20635500.0    4659000.00
pj_inv result: 118.622841      40.452777

 

 注：由于我国北京54和西安80中的地理坐标系的基准数据没有公开，因此无法完成准确计算。