第11章-坐标系转换与投影

第11章：坐标系转换与投影

11.1 坐标系统概述

坐标系统是 GIS 的基础，用于精确描述地球上位置的数学模型。不同的应用场景需要使用不同的坐标系统。

11.1.1 坐标系统类型

类型	说明	示例
地理坐标系	基于椭球体，使用经纬度	WGS84 (EPSG:4326)
投影坐标系	将地球投影到平面，使用米/英尺	UTM, Web Mercator
本地坐标系	局部平面坐标系	工程测量坐标

11.1.2 常用坐标系统

EPSG 代码	名称	说明
4326	WGS84	GPS 标准坐标系
3857	Web Mercator	Web 地图标准
4490	CGCS2000	中国大地坐标系
4547	CGCS2000 / 3-degree Gauss-Kruger zone 38	中国高斯投影
32650	WGS 84 / UTM zone 50N	UTM 50北带

11.1.3 安装 ProjNet

dotnet add package ProjNet

11.2 ProjNet 基础

11.2.1 创建坐标系统

using ProjNet.CoordinateSystems;

var csFactory = new CoordinateSystemFactory();

// 从 WKT 创建坐标系统
var wgs84Wkt = @"
    GEOGCS[""WGS 84"",
        DATUM[""WGS_1984"",
            SPHEROID[""WGS 84"",6378137,298.257223563]],
        PRIMEM[""Greenwich"",0],
        UNIT[""degree"",0.0174532925199433]]";

var wgs84 = csFactory.CreateFromWkt(wgs84Wkt) as GeographicCoordinateSystem;
Console.WriteLine($"坐标系名称: {wgs84.Name}");
Console.WriteLine($"椭球体: {wgs84.HorizontalDatum.Ellipsoid.Name}");

11.2.2 预定义坐标系统

using ProjNet.CoordinateSystems;

// 创建常用坐标系统
var gcsFactory = new GeographicCoordinateSystem(
    "WGS 84",
    AngularUnit.Degrees,
    HorizontalDatum.WGS84,
    PrimeMeridian.Greenwich,
    new AxisInfo("Lon", AxisOrientationEnum.East),
    new AxisInfo("Lat", AxisOrientationEnum.North));

// Web Mercator
var webMercatorWkt = @"
    PROJCS[""WGS 84 / Pseudo-Mercator"",
        GEOGCS[""WGS 84"",
            DATUM[""WGS_1984"",
                SPHEROID[""WGS 84"",6378137,298.257223563]],
            PRIMEM[""Greenwich"",0],
            UNIT[""degree"",0.0174532925199433]],
        PROJECTION[""Mercator_1SP""],
        PARAMETER[""central_meridian"",0],
        PARAMETER[""scale_factor"",1],
        PARAMETER[""false_easting"",0],
        PARAMETER[""false_northing"",0],
        UNIT[""metre"",1]]";

var webMercator = csFactory.CreateFromWkt(webMercatorWkt) as ProjectedCoordinateSystem;

11.2.3 使用 SRID 服务

using ProjNet.CoordinateSystems;
using ProjNet.IO.CoordinateSystems;

// 创建 SRID 字典
var sridReader = new SridReader();

// 从文件加载
// sridReader.Load("srid.csv");

// 或从资源加载常用坐标系
public static class CommonCoordinateSystems
{
    private static readonly CoordinateSystemFactory Factory = new();
    
    public static CoordinateSystem WGS84 => Factory.CreateFromWkt(@"
        GEOGCS[""WGS 84"",
            DATUM[""WGS_1984"",SPHEROID[""WGS 84"",6378137,298.257223563]],
            PRIMEM[""Greenwich"",0],
            UNIT[""degree"",0.0174532925199433]]");
    
    public static CoordinateSystem WebMercator => Factory.CreateFromWkt(@"
        PROJCS[""WGS 84 / Pseudo-Mercator"",
            GEOGCS[""WGS 84"",
                DATUM[""WGS_1984"",SPHEROID[""WGS 84"",6378137,298.257223563]],
                PRIMEM[""Greenwich"",0],
                UNIT[""degree"",0.0174532925199433]],
            PROJECTION[""Mercator_1SP""],
            PARAMETER[""central_meridian"",0],
            PARAMETER[""scale_factor"",1],
            PARAMETER[""false_easting"",0],
            PARAMETER[""false_northing"",0],
            UNIT[""metre"",1]]");
    
    public static CoordinateSystem CGCS2000 => Factory.CreateFromWkt(@"
        GEOGCS[""China Geodetic Coordinate System 2000"",
            DATUM[""China_2000"",
                SPHEROID[""CGCS2000"",6378137,298.257222101]],
            PRIMEM[""Greenwich"",0],
            UNIT[""degree"",0.0174532925199433]]");
}

11.3 坐标转换

11.3.1 基本坐标转换

using ProjNet.CoordinateSystems;
using ProjNet.CoordinateSystems.Transformations;

var csFactory = new CoordinateSystemFactory();
var ctFactory = new CoordinateTransformationFactory();

// 创建源坐标系（WGS84）
var wgs84 = csFactory.CreateFromWkt(@"
    GEOGCS[""WGS 84"",
        DATUM[""WGS_1984"",SPHEROID[""WGS 84"",6378137,298.257223563]],
        PRIMEM[""Greenwich"",0],
        UNIT[""degree"",0.0174532925199433]]");

// 创建目标坐标系（UTM Zone 50N）
var utm50n = csFactory.CreateFromWkt(@"
    PROJCS[""WGS 84 / UTM zone 50N"",
        GEOGCS[""WGS 84"",
            DATUM[""WGS_1984"",SPHEROID[""WGS 84"",6378137,298.257223563]],
            PRIMEM[""Greenwich"",0],
            UNIT[""degree"",0.0174532925199433]],
        PROJECTION[""Transverse_Mercator""],
        PARAMETER[""latitude_of_origin"",0],
        PARAMETER[""central_meridian"",117],
        PARAMETER[""scale_factor"",0.9996],
        PARAMETER[""false_easting"",500000],
        PARAMETER[""false_northing"",0],
        UNIT[""metre"",1]]");

// 创建转换
var transform = ctFactory.CreateFromCoordinateSystems(wgs84, utm50n);

// 转换坐标（北京）
var wgs84Coord = new double[] { 116.4074, 39.9042 };
var utmCoord = transform.MathTransform.Transform(wgs84Coord);

Console.WriteLine($"WGS84: ({wgs84Coord[0]}, {wgs84Coord[1]})");
Console.WriteLine($"UTM: ({utmCoord[0]:F2}, {utmCoord[1]:F2})");

11.3.2 批量坐标转换

// 批量转换
var coordinates = new double[][]
{
    new double[] { 116.4074, 39.9042 },  // 北京
    new double[] { 121.4737, 31.2304 },  // 上海
    new double[] { 113.2644, 23.1291 }   // 广州
};

var transformedCoords = new List<double[]>();
foreach (var coord in coordinates)
{
    var transformed = transform.MathTransform.Transform(coord);
    transformedCoords.Add(transformed);
    Console.WriteLine($"({coord[0]}, {coord[1]}) -> ({transformed[0]:F2}, {transformed[1]:F2})");
}

// 逆向转换
var inverseTransform = transform.MathTransform.Inverse();
foreach (var coord in transformedCoords)
{
    var original = inverseTransform.Transform(coord);
    Console.WriteLine($"逆转换: ({coord[0]:F2}, {coord[1]:F2}) -> ({original[0]:F6}, {original[1]:F6})");
}

11.3.3 与 NTS 几何集成

using NetTopologySuite.Geometries;
using NetTopologySuite.Geometries.Utilities;
using ProjNet.CoordinateSystems.Transformations;

public class CoordinateTransformFilter : ICoordinateSequenceFilter
{
    private readonly MathTransform _transform;
    
    public bool Done => false;
    public bool GeometryChanged => true;
    
    public CoordinateTransformFilter(MathTransform transform)
    {
        _transform = transform;
    }
    
    public void Filter(CoordinateSequence seq, int i)
    {
        var x = seq.GetX(i);
        var y = seq.GetY(i);
        var transformed = _transform.Transform(new double[] { x, y });
        seq.SetX(i, transformed[0]);
        seq.SetY(i, transformed[1]);
    }
}

// 使用示例
var factory = new GeometryFactory(new PrecisionModel(), 4326);

var polygon = factory.CreatePolygon(new Coordinate[]
{
    new Coordinate(116.3, 39.8), new Coordinate(116.5, 39.8),
    new Coordinate(116.5, 40.0), new Coordinate(116.3, 40.0),
    new Coordinate(116.3, 39.8)
});

Console.WriteLine($"转换前: {polygon.AsText()}");

// 创建转换
var transform = ctFactory.CreateFromCoordinateSystems(wgs84, utm50n);

// 复制几何并转换
var transformedPolygon = (Polygon)polygon.Copy();
transformedPolygon.Apply(new CoordinateTransformFilter(transform.MathTransform));

// 更新 SRID
transformedPolygon = new GeometryFactory(new PrecisionModel(), 32650)
    .CreatePolygon(transformedPolygon.Coordinates);

Console.WriteLine($"转换后: {transformedPolygon.AsText()}");
Console.WriteLine($"SRID: {transformedPolygon.SRID}");

11.4 坐标转换服务

11.4.1 完整的转换服务

using NetTopologySuite.Geometries;
using ProjNet.CoordinateSystems;
using ProjNet.CoordinateSystems.Transformations;

public class CoordinateTransformService
{
    private readonly CoordinateSystemFactory _csFactory;
    private readonly CoordinateTransformationFactory _ctFactory;
    private readonly Dictionary<int, CoordinateSystem> _coordinateSystems;

    public CoordinateTransformService()
    {
        _csFactory = new CoordinateSystemFactory();
        _ctFactory = new CoordinateTransformationFactory();
        _coordinateSystems = new Dictionary<int, CoordinateSystem>();
        
        InitializeCommonCoordinateSystems();
    }

    private void InitializeCommonCoordinateSystems()
    {
        // WGS84
        _coordinateSystems[4326] = _csFactory.CreateFromWkt(@"
            GEOGCS[""WGS 84"",
                DATUM[""WGS_1984"",SPHEROID[""WGS 84"",6378137,298.257223563]],
                PRIMEM[""Greenwich"",0],
                UNIT[""degree"",0.0174532925199433]]");

        // Web Mercator
        _coordinateSystems[3857] = _csFactory.CreateFromWkt(@"
            PROJCS[""WGS 84 / Pseudo-Mercator"",
                GEOGCS[""WGS 84"",
                    DATUM[""WGS_1984"",SPHEROID[""WGS 84"",6378137,298.257223563]],
                    PRIMEM[""Greenwich"",0],
                    UNIT[""degree"",0.0174532925199433]],
                PROJECTION[""Mercator_1SP""],
                PARAMETER[""central_meridian"",0],
                PARAMETER[""scale_factor"",1],
                PARAMETER[""false_easting"",0],
                PARAMETER[""false_northing"",0],
                UNIT[""metre"",1]]");

        // CGCS2000
        _coordinateSystems[4490] = _csFactory.CreateFromWkt(@"
            GEOGCS[""China Geodetic Coordinate System 2000"",
                DATUM[""China_2000"",SPHEROID[""CGCS2000"",6378137,298.257222101]],
                PRIMEM[""Greenwich"",0],
                UNIT[""degree"",0.0174532925199433]]");
    }

    /// <summary>
    /// 注册坐标系统
    /// </summary>
    public void RegisterCoordinateSystem(int srid, string wkt)
    {
        _coordinateSystems[srid] = _csFactory.CreateFromWkt(wkt);
    }

    /// <summary>
    /// 转换坐标点
    /// </summary>
    public (double X, double Y) TransformPoint(
        double x, double y, 
        int sourceSrid, int targetSrid)
    {
        if (sourceSrid == targetSrid)
            return (x, y);

        var sourceCs = _coordinateSystems[sourceSrid];
        var targetCs = _coordinateSystems[targetSrid];
        
        var transform = _ctFactory.CreateFromCoordinateSystems(sourceCs, targetCs);
        var result = transform.MathTransform.Transform(new double[] { x, y });
        
        return (result[0], result[1]);
    }

    /// <summary>
    /// 转换几何对象
    /// </summary>
    public Geometry TransformGeometry(Geometry geometry, int targetSrid)
    {
        if (geometry.SRID == targetSrid)
            return geometry;

        var sourceCs = _coordinateSystems[geometry.SRID];
        var targetCs = _coordinateSystems[targetSrid];
        
        var transform = _ctFactory.CreateFromCoordinateSystems(sourceCs, targetCs);
        
        // 复制几何
        var transformed = (Geometry)geometry.Copy();
        
        // 应用转换
        transformed.Apply(new CoordinateTransformFilter(transform.MathTransform));
        
        // 创建新几何（带有正确的 SRID）
        var targetFactory = new GeometryFactory(
            new PrecisionModel(), targetSrid);
        
        return targetFactory.CreateGeometry(transformed);
    }

    /// <summary>
    /// 批量转换几何
    /// </summary>
    public IEnumerable<Geometry> TransformGeometries(
        IEnumerable<Geometry> geometries, 
        int targetSrid)
    {
        return geometries.Select(g => TransformGeometry(g, targetSrid));
    }

    /// <summary>
    /// 计算真实距离（米）
    /// </summary>
    public double CalculateDistance(
        double lon1, double lat1, 
        double lon2, double lat2)
    {
        // 使用 Haversine 公式
        const double R = 6371000; // 地球半径（米）
        
        var lat1Rad = lat1 * Math.PI / 180;
        var lat2Rad = lat2 * Math.PI / 180;
        var deltaLat = (lat2 - lat1) * Math.PI / 180;
        var deltaLon = (lon2 - lon1) * Math.PI / 180;
        
        var a = Math.Sin(deltaLat / 2) * Math.Sin(deltaLat / 2) +
                Math.Cos(lat1Rad) * Math.Cos(lat2Rad) *
                Math.Sin(deltaLon / 2) * Math.Sin(deltaLon / 2);
        
        var c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
        
        return R * c;
    }

    /// <summary>
    /// WGS84 转 Web Mercator
    /// </summary>
    public (double X, double Y) Wgs84ToWebMercator(double lon, double lat)
    {
        return TransformPoint(lon, lat, 4326, 3857);
    }

    /// <summary>
    /// Web Mercator 转 WGS84
    /// </summary>
    public (double Lon, double Lat) WebMercatorToWgs84(double x, double y)
    {
        return TransformPoint(x, y, 3857, 4326);
    }
}

11.4.2 使用示例

var transformService = new CoordinateTransformService();

// 坐标转换
var (x, y) = transformService.Wgs84ToWebMercator(116.4074, 39.9042);
Console.WriteLine($"北京 Web Mercator: ({x:F2}, {y:F2})");

var (lon, lat) = transformService.WebMercatorToWgs84(x, y);
Console.WriteLine($"转回 WGS84: ({lon:F6}, {lat:F6})");

// 几何转换
var factory = new GeometryFactory(new PrecisionModel(), 4326);
var polygon = factory.CreatePolygon(new Coordinate[]
{
    new Coordinate(116.3, 39.8), new Coordinate(116.5, 39.8),
    new Coordinate(116.5, 40.0), new Coordinate(116.3, 40.0),
    new Coordinate(116.3, 39.8)
});

var projectedPolygon = transformService.TransformGeometry(polygon, 3857);
Console.WriteLine($"投影后面积: {projectedPolygon.Area:F2} 平方米");

// 计算真实距离
var distance = transformService.CalculateDistance(
    116.4074, 39.9042,  // 北京
    121.4737, 31.2304   // 上海
);
Console.WriteLine($"北京到上海距离: {distance / 1000:F2} 公里");

11.5 中国常用坐标系转换

11.5.1 WGS84、GCJ02、BD09 互转

/// <summary>
/// 中国常用坐标系转换（WGS84、GCJ02、BD09）
/// </summary>
public static class ChineseCoordinateTransform
{
    private const double PI = Math.PI;
    private const double X_PI = PI * 3000.0 / 180.0;
    private const double A = 6378245.0;
    private const double EE = 0.00669342162296594323;

    /// <summary>
    /// WGS84 转 GCJ02（火星坐标）
    /// </summary>
    public static (double Lon, double Lat) Wgs84ToGcj02(double lon, double lat)
    {
        if (OutOfChina(lon, lat))
            return (lon, lat);

        var dLat = TransformLat(lon - 105.0, lat - 35.0);
        var dLon = TransformLon(lon - 105.0, lat - 35.0);
        var radLat = lat / 180.0 * PI;
        var magic = Math.Sin(radLat);
        magic = 1 - EE * magic * magic;
        var sqrtMagic = Math.Sqrt(magic);
        dLat = (dLat * 180.0) / ((A * (1 - EE)) / (magic * sqrtMagic) * PI);
        dLon = (dLon * 180.0) / (A / sqrtMagic * Math.Cos(radLat) * PI);

        return (lon + dLon, lat + dLat);
    }

    /// <summary>
    /// GCJ02 转 WGS84
    /// </summary>
    public static (double Lon, double Lat) Gcj02ToWgs84(double lon, double lat)
    {
        if (OutOfChina(lon, lat))
            return (lon, lat);

        var (mLon, mLat) = Wgs84ToGcj02(lon, lat);
        return (lon * 2 - mLon, lat * 2 - mLat);
    }

    /// <summary>
    /// GCJ02 转 BD09（百度坐标）
    /// </summary>
    public static (double Lon, double Lat) Gcj02ToBd09(double lon, double lat)
    {
        var z = Math.Sqrt(lon * lon + lat * lat) + 0.00002 * Math.Sin(lat * X_PI);
        var theta = Math.Atan2(lat, lon) + 0.000003 * Math.Cos(lon * X_PI);
        return (z * Math.Cos(theta) + 0.0065, z * Math.Sin(theta) + 0.006);
    }

    /// <summary>
    /// BD09 转 GCJ02
    /// </summary>
    public static (double Lon, double Lat) Bd09ToGcj02(double lon, double lat)
    {
        var x = lon - 0.0065;
        var y = lat - 0.006;
        var z = Math.Sqrt(x * x + y * y) - 0.00002 * Math.Sin(y * X_PI);
        var theta = Math.Atan2(y, x) - 0.000003 * Math.Cos(x * X_PI);
        return (z * Math.Cos(theta), z * Math.Sin(theta));
    }

    /// <summary>
    /// WGS84 转 BD09
    /// </summary>
    public static (double Lon, double Lat) Wgs84ToBd09(double lon, double lat)
    {
        var gcj = Wgs84ToGcj02(lon, lat);
        return Gcj02ToBd09(gcj.Lon, gcj.Lat);
    }

    /// <summary>
    /// BD09 转 WGS84
    /// </summary>
    public static (double Lon, double Lat) Bd09ToWgs84(double lon, double lat)
    {
        var gcj = Bd09ToGcj02(lon, lat);
        return Gcj02ToWgs84(gcj.Lon, gcj.Lat);
    }

    private static bool OutOfChina(double lon, double lat)
    {
        return lon < 72.004 || lon > 137.8347 || lat < 0.8293 || lat > 55.8271;
    }

    private static double TransformLat(double x, double y)
    {
        var ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 
                  0.2 * Math.Sqrt(Math.Abs(x));
        ret += (20.0 * Math.Sin(6.0 * x * PI) + 20.0 * Math.Sin(2.0 * x * PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.Sin(y * PI) + 40.0 * Math.Sin(y / 3.0 * PI)) * 2.0 / 3.0;
        ret += (160.0 * Math.Sin(y / 12.0 * PI) + 320 * Math.Sin(y * PI / 30.0)) * 2.0 / 3.0;
        return ret;
    }

    private static double TransformLon(double x, double y)
    {
        var ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 
                  0.1 * Math.Sqrt(Math.Abs(x));
        ret += (20.0 * Math.Sin(6.0 * x * PI) + 20.0 * Math.Sin(2.0 * x * PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.Sin(x * PI) + 40.0 * Math.Sin(x / 3.0 * PI)) * 2.0 / 3.0;
        ret += (150.0 * Math.Sin(x / 12.0 * PI) + 300.0 * Math.Sin(x / 30.0 * PI)) * 2.0 / 3.0;
        return ret;
    }
}

// 使用示例
var (gcjLon, gcjLat) = ChineseCoordinateTransform.Wgs84ToGcj02(116.4074, 39.9042);
Console.WriteLine($"GCJ02: ({gcjLon:F6}, {gcjLat:F6})");

var (bdLon, bdLat) = ChineseCoordinateTransform.Wgs84ToBd09(116.4074, 39.9042);
Console.WriteLine($"BD09: ({bdLon:F6}, {bdLat:F6})");

11.6 本章小结

本章详细介绍了坐标系转换与投影：

1. 坐标系统基础：地理坐标系、投影坐标系
2. ProjNet 使用：创建坐标系统、坐标转换
3. 与 NTS 集成：几何对象的坐标转换
4. 转换服务：完整的坐标转换服务类
5. 中国坐标系：WGS84、GCJ02、BD09 互转

11.7 下一步

下一章我们将学习矢量切片生成。

---

相关资源：

• ProjNet GitHub
• EPSG.io - 坐标系统查询
• spatialreference.org