11.6 本章小结

副本挑战

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

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

public CoordinateTransformService()

{

_csFactory = new CoordinateSystemFactory();

_ctFactory = new CoordinateTransformationFactory();

_coordinateSystems = new Dictionary();

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

}

///

/// 注册坐标系统

///

public void RegisterCoordinateSystem(int srid, string wkt)

{

_coordinateSystems[srid] = _csFactory.CreateFromWkt(wkt);

}

///

/// 转换坐标点

///

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

}

///

/// 转换几何对象

///

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

}

///

/// 批量转换几何

///

public IEnumerable TransformGeometries(

IEnumerable geometries,

int targetSrid)

{

return geometries.Select(g => TransformGeometry(g, targetSrid));

}

///

/// 计算真实距离（米）

///

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;

}

///

/// WGS84 转 Web Mercator

///

public (double X, double Y) Wgs84ToWebMercator(double lon, double lat)

{

return TransformPoint(lon, lat, 4326, 3857);

}

///

/// Web Mercator 转 WGS84

///

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 互转

///

/// 中国常用坐标系转换（WGS84、GCJ02、BD09）

///

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;

///

/// WGS84 转 GCJ02（火星坐标）

///

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

}

///

/// GCJ02 转 WGS84

///

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

}

///

/// GCJ02 转 BD09（百度坐标）

///

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

}

///

/// BD09 转 GCJ02

///

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

}

///

/// WGS84 转 BD09

///

public static (double Lon, double Lat) Wgs84ToBd09(double lon, double lat)

{

var gcj = Wgs84ToGcj02(lon, lat);

return Gcj02ToBd09(gcj.Lon, gcj.Lat);

}

///

/// BD09 转 WGS84

///

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 本章小结

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

坐标系统基础：地理坐标系、投影坐标系

ProjNet 使用：创建坐标系统、坐标转换

与 NTS 集成：几何对象的坐标转换

转换服务：完整的坐标转换服务类

中国坐标系：WGS84、GCJ02、BD09 互转

11.7 下一步

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

相关资源：

ProjNet GitHub

EPSG.io - 坐标系统查询

spatialreference.org