第08章 - 坐标系管理与转换
第08章 - 坐标系管理与转换
8.1 坐标系基础
8.1.1 坐标系分类
GIS中常用的坐标系分为两大类：
地理坐标系 (Geographic Coordinate System, GCS)
- 使用经度和纬度表示位置
- 单位是度
- 常见：WGS84 (EPSG:4326), CGCS2000 (EPSG:4490)
投影坐标系 (Projected Coordinate System, PCS)
- 将球面坐标投影到平面
- 单位通常是米
- 常见：UTM, 高斯-克吕格投影
8.1.2 常用坐标系
| 坐标系 | EPSG代码 | 类型 | 说明 |
|---|---|---|---|
| WGS 84 | 4326 | 地理 | GPS全球定位系统 |
| CGCS2000 | 4490 | 地理 | 中国大地坐标系2000 |
| WGS 84 / UTM zone 50N | 32650 | 投影 | UTM 50带北半球 |
| CGCS2000 / 3-degree Zone 39 | 4520 | 投影 | CGCS2000 3度带39带 |
8.1.3 中国坐标系特点
中国主要使用CGCS2000（2000国家大地坐标系）：
- CGCS2000地理坐标系：EPSG:4490
- CGCS2000投影坐标系：
- 3度带：EPSG:4491-4554（带号24-45）
- 6度带：EPSG:4513-4533（带号13-23）
8.2 CrsUtil工具类
8.2.1 功能概述
CrsUtil 是OGU4Net的坐标系工具类，提供：
| 功能 | 方法 | 说明 |
|---|---|---|
| 坐标转换 | Transform | WKT或Geometry坐标转换 |
| 带号计算 | GetDh, GetDh6 | 根据经度计算带号 |
| WKID转换 | GetProjectedWkid | 根据带号获取投影WKID |
| 容差获取 | GetTolerance | 获取坐标系适合的容差 |
| 类型判断 | IsProjectedCRS | 判断是否为投影坐标系 |
8.2.2 基本使用
using OpenGIS.Utils.Engine.Util;
// 坐标转换
string wkt = "POINT (116.404 39.915)";
string transformed = CrsUtil.Transform(wkt, 4326, 4490);
Console.WriteLine($"WGS84 → CGCS2000: {transformed}");
// 获取带号
double longitude = 116.404;
int zone3 = CrsUtil.GetDh(longitude); // 3度带带号
int zone6 = CrsUtil.GetDh6(longitude); // 6度带带号
Console.WriteLine($"经度 {longitude}: 3度带第{zone3}带, 6度带第{zone6}带");
8.3 坐标转换
8.3.1 WKT坐标转换
// WGS84 (4326) → CGCS2000 (4490)
string wgs84Point = "POINT (116.404 39.915)";
string cgcs2000Point = CrsUtil.Transform(wgs84Point, 4326, 4490);
Console.WriteLine(cgcs2000Point);
// WGS84 → CGCS2000投影坐标系（3度带39带）
string projected = CrsUtil.Transform(wgs84Point, 4326, 4520);
Console.WriteLine(projected); // 坐标变为米单位
// 多边形转换
string wgs84Polygon = "POLYGON ((116.3 39.8, 116.5 39.8, 116.5 40.0, 116.3 40.0, 116.3 39.8))";
string cgcs2000Polygon = CrsUtil.Transform(wgs84Polygon, 4326, 4490);
8.3.2 Geometry对象转换
using OpenGIS.Utils.Geometry;
using OgrGeometry = OSGeo.OGR.Geometry;
// 创建Geometry对象
OgrGeometry geom = GeometryUtil.Wkt2Geometry("POINT (116.404 39.915)");
// 转换坐标系
OgrGeometry transformed = CrsUtil.Transform(geom, 4326, 4490);
// 转回WKT
string wkt = GeometryUtil.Geometry2Wkt(transformed);
Console.WriteLine(wkt);
8.3.3 批量坐标转换
// 转换整个图层的坐标
var layer = OguLayerUtil.ReadLayer(DataFormatType.SHP, "wgs84.shp");
foreach (var feature in layer.Features)
{
if (!string.IsNullOrEmpty(feature.Wkt))
{
feature.Wkt = CrsUtil.Transform(feature.Wkt, 4326, 4490);
}
}
// 更新图层坐标系信息
layer.Wkid = 4490;
// 保存
OguLayerUtil.WriteLayer(DataFormatType.SHP, layer, "cgcs2000.shp");
8.3.4 带投影的转换
// 经纬度转投影坐标（适合面积计算）
string geoWkt = "POLYGON ((116.3 39.8, 116.5 39.8, 116.5 40.0, 116.3 40.0, 116.3 39.8))";
// 计算带号
int zone = CrsUtil.GetDh(116.4); // 39带
// 获取对应的投影WKID
int projWkid = CrsUtil.GetProjectedWkid(zone); // 4520
// 转换到投影坐标系
string projWkt = CrsUtil.Transform(geoWkt, 4326, projWkid);
// 计算面积（单位：平方米）
double area = GeometryUtil.AreaWkt(projWkt);
Console.WriteLine($"面积: {area / 10000:F2} 公顷");
8.4 带号计算
8.4.1 3度带计算
// 中国3度带范围：24带至45带
// 中央经线 = 带号 × 3
// 根据经度计算3度带带号
double[] longitudes = { 73.0, 90.0, 105.0, 116.4, 120.0, 135.0 };
foreach (var lon in longitudes)
{
int zone = CrsUtil.GetDh(lon);
int centralMeridian = zone * 3;
Console.WriteLine($"经度 {lon}° → 3度带第{zone}带 (中央经线 {centralMeridian}°)");
}
// 输出示例：
// 经度 116.4° → 3度带第39带 (中央经线 117°)
8.4.2 6度带计算
// 中国6度带范围：13带至23带
// 中央经线 = 带号 × 6 - 3
double[] longitudes = { 75.0, 90.0, 105.0, 116.4, 120.0, 135.0 };
foreach (var lon in longitudes)
{
int zone = CrsUtil.GetDh6(lon);
int centralMeridian = zone * 6 - 3;
Console.WriteLine($"经度 {lon}° → 6度带第{zone}带 (中央经线 {centralMeridian}°)");
}
// 输出示例：
// 经度 116.4° → 6度带第20带 (中央经线 117°)
8.4.3 根据Geometry计算带号
// 根据几何对象的质心计算带号
var geom = GeometryUtil.Wkt2Geometry(
"POLYGON ((116.3 39.8, 116.5 39.8, 116.5 40.0, 116.3 40.0, 116.3 39.8))");
int zone = CrsUtil.GetDh(geom);
Console.WriteLine($"几何对象所在带号: {zone}");
8.5 WKID转换
8.5.1 带号转WKID
// 3度带：带号 → WKID
// CGCS2000 3度带WKID范围：4491-4554
int zone3 = 39;
int wkid3 = CrsUtil.GetProjectedWkid(zone3);
Console.WriteLine($"3度带第{zone3}带 → EPSG:{wkid3}"); // EPSG:4520
// 6度带：带号 → WKID
// CGCS2000 6度带WKID范围：4513-4533
int zone6 = 20;
int wkid6 = CrsUtil.GetProjectedWkid6(zone6);
Console.WriteLine($"6度带第{zone6}带 → EPSG:{wkid6}"); // EPSG:4520
8.5.2 WKID转带号
// WKID → 带号
int wkid = 4520;
int zone = CrsUtil.GetDhFromWkid(wkid);
Console.WriteLine($"EPSG:{wkid} → 第{zone}带");
8.5.3 中国坐标系WKID对照表
CGCS2000 3度带（常用范围）：
| 带号 | 中央经线 | EPSG代码 | 适用区域 |
|---|---|---|---|
| 35 | 105° | 4516 | 四川、云南东部 |
| 36 | 108° | 4517 | 陕西、重庆 |
| 37 | 111° | 4518 | 河南、湖北 |
| 38 | 114° | 4519 | 河北、山东西部 |
| 39 | 117° | 4520 | 北京、天津、山东东部 |
| 40 | 120° | 4521 | 上海、江苏、浙江 |
| 41 | 123° | 4522 | 辽宁、吉林南部 |
8.6 容差与精度
8.6.1 获取推荐容差
// 不同坐标系需要不同的容差值
// 地理坐标系（度）
double tolGeo = CrsUtil.GetTolerance(4326);
Console.WriteLine($"WGS84容差: {tolGeo}"); // 0.0000001
double tolCgcs = CrsUtil.GetTolerance(4490);
Console.WriteLine($"CGCS2000容差: {tolCgcs}"); // 0.0000001
// 投影坐标系（米）
double tolProj = CrsUtil.GetTolerance(4520);
Console.WriteLine($"投影坐标系容差: {tolProj}"); // 0.001
8.6.2 坐标精度处理
// 经纬度坐标一般保留6-8位小数
// 6位小数 ≈ 0.1米精度
// 8位小数 ≈ 1毫米精度
double lat = 39.9150678901234;
double lon = 116.4039876543210;
// 保留6位小数
lat = Math.Round(lat, 6); // 39.915068
lon = Math.Round(lon, 6); // 116.403988
8.7 判断坐标系类型
8.7.1 判断是否为投影坐标系
// 判断WKID是否为投影坐标系
bool isProj4326 = CrsUtil.IsProjectedCRS(4326);
Console.WriteLine($"4326是投影坐标系: {isProj4326}"); // false
bool isProj4490 = CrsUtil.IsProjectedCRS(4490);
Console.WriteLine($"4490是投影坐标系: {isProj4490}"); // false
bool isProj4520 = CrsUtil.IsProjectedCRS(4520);
Console.WriteLine($"4520是投影坐标系: {isProj4520}"); // true
bool isProj32650 = CrsUtil.IsProjectedCRS(32650);
Console.WriteLine($"32650是投影坐标系: {isProj32650}"); // true
8.7.2 根据坐标系选择处理方式
public double CalculateArea(OguLayer layer)
{
double totalArea = 0;
foreach (var feature in layer.Features)
{
if (string.IsNullOrEmpty(feature.Wkt)) continue;
string wkt = feature.Wkt;
// 如果是地理坐标系，需要先转换到投影坐标系
if (layer.Wkid.HasValue && !CrsUtil.IsProjectedCRS(layer.Wkid.Value))
{
// 计算带号
var geom = GeometryUtil.Wkt2Geometry(wkt);
int zone = CrsUtil.GetDh(geom);
int projWkid = CrsUtil.GetProjectedWkid(zone);
// 转换
wkt = CrsUtil.Transform(wkt, layer.Wkid.Value, projWkid);
}
totalArea += GeometryUtil.AreaWkt(wkt);
}
return totalArea;
}
8.8 实战应用
8.8.1 GPS数据转换
// GPS设备输出的是WGS84坐标
var gpsPoints = new[]
{
("POINT (116.404 39.915)", "点1"),
("POINT (116.405 39.916)", "点2"),
("POINT (116.406 39.917)", "点3")
};
// 创建CGCS2000图层
var layer = new OguLayer
{
Name = "GPS测量点",
GeometryType = GeometryType.POINT,
Wkid = 4490
};
layer.AddField(new OguField { Name = "Name", DataType = FieldDataType.STRING, Length = 50 });
int fid = 1;
foreach (var (wkt, name) in gpsPoints)
{
// WGS84 → CGCS2000
var cgcsWkt = CrsUtil.Transform(wkt, 4326, 4490);
var feature = new OguFeature { Fid = fid++, Wkt = cgcsWkt };
feature.SetValue("Name", name);
layer.AddFeature(feature);
}
OguLayerUtil.WriteLayer(DataFormatType.SHP, layer, "gps_cgcs2000.shp");
8.8.2 面积计算标准化
// 计算地块面积（结果单位：公顷）
public double CalculatePlotArea(string wkt, int wkid)
{
string projWkt;
if (CrsUtil.IsProjectedCRS(wkid))
{
// 已经是投影坐标系
projWkt = wkt;
}
else
{
// 地理坐标系，需要转换
var geom = GeometryUtil.Wkt2Geometry(wkt);
int zone = CrsUtil.GetDh(geom);
int projWkid = CrsUtil.GetProjectedWkid(zone);
projWkt = CrsUtil.Transform(wkt, wkid, projWkid);
}
double areaM2 = GeometryUtil.AreaWkt(projWkt);
return areaM2 / 10000; // 转换为公顷
}
// 使用示例
var wkt = "POLYGON ((116.3 39.8, 116.5 39.8, 116.5 40.0, 116.3 40.0, 116.3 39.8))";
double area = CalculatePlotArea(wkt, 4326);
Console.WriteLine($"面积: {area:F2} 公顷");
8.8.3 跨带处理
// 处理跨越两个带的数据
public OguLayer ProcessCrossZoneLayer(OguLayer layer)
{
// 按带号分组
var featuresByZone = layer.Features
.GroupBy(f =>
{
if (string.IsNullOrEmpty(f.Wkt)) return 0;
var geom = GeometryUtil.Wkt2Geometry(f.Wkt);
return CrsUtil.GetDh(geom);
})
.Where(g => g.Key > 0)
.ToDictionary(g => g.Key, g => g.ToList());
Console.WriteLine($"数据跨越 {featuresByZone.Count} 个带:");
foreach (var (zone, features) in featuresByZone)
{
Console.WriteLine($" 第{zone}带: {features.Count} 个要素");
}
// 如果跨带，使用最多要素的带号
int mainZone = featuresByZone
.OrderByDescending(kvp => kvp.Value.Count)
.First().Key;
int projWkid = CrsUtil.GetProjectedWkid(mainZone);
Console.WriteLine($"使用第{mainZone}带 (EPSG:{projWkid}) 作为统一投影");
// 转换所有要素
var result = layer.Clone();
result.Wkid = projWkid;
foreach (var feature in result.Features)
{
if (!string.IsNullOrEmpty(feature.Wkt))
{
feature.Wkt = CrsUtil.Transform(feature.Wkt, layer.Wkid ?? 4490, projWkid);
}
}
return result;
}
8.8.4 国土数据转换
// 国土部门经常使用CGCS2000投影坐标系
// 需要转换为地理坐标系用于Web展示
var layer = OguLayerUtil.ReadLayer(DataFormatType.SHP, "cadastre.shp");
// 假设原始数据是CGCS2000 3度带39带
int srcWkid = 4520;
// 转换为CGCS2000地理坐标系
foreach (var feature in layer.Features)
{
if (!string.IsNullOrEmpty(feature.Wkt))
{
feature.Wkt = CrsUtil.Transform(feature.Wkt, srcWkid, 4490);
}
}
layer.Wkid = 4490;
// 输出为GeoJSON供Web使用
// Web地图通常需要WGS84，再转换一次
foreach (var feature in layer.Features)
{
if (!string.IsNullOrEmpty(feature.Wkt))
{
feature.Wkt = CrsUtil.Transform(feature.Wkt, 4490, 4326);
}
}
layer.Wkid = 4326;
OguLayerUtil.WriteLayer(DataFormatType.GEOJSON, layer, "cadastre_web.geojson");
8.9 注意事项
8.9.1 精度损失
// 多次转换会累积精度损失
// 尽量减少转换次数
// 不推荐：多次转换
var wkt1 = CrsUtil.Transform(wkt, 4326, 4490);
var wkt2 = CrsUtil.Transform(wkt1, 4490, 4520);
var wkt3 = CrsUtil.Transform(wkt2, 4520, 4326);
// 推荐：一步到位
var wktDirect = CrsUtil.Transform(wkt, 4326, 4520);
8.9.2 大数据集性能
// 对于大数据集，考虑使用GDAL命令行工具
// ogr2ogr -t_srs EPSG:4490 output.shp input.shp
// 或者使用并行处理
var features = layer.Features.ToList();
Parallel.ForEach(features, feature =>
{
if (!string.IsNullOrEmpty(feature.Wkt))
{
feature.Wkt = CrsUtil.Transform(feature.Wkt, 4326, 4490);
}
});
8.10 小结
本章介绍了OGU4Net的坐标系管理功能：
- 坐标系基础：地理坐标系和投影坐标系的区别
- CrsUtil工具类：坐标转换、带号计算、WKID转换
- 坐标转换：WKT和Geometry对象的坐标系转换
- 带号计算：3度带和6度带的计算方法
- WKID转换：带号与EPSG代码的相互转换
- 实战应用：GPS数据转换、面积计算、跨带处理
掌握坐标系转换是GIS开发的基础技能，正确处理坐标系可以避免很多数据问题。
浙公网安备 33010602011771号