坐标参考系统与投影

第八章 坐标参考系统与投影

8.1 坐标参考系统基础

8.1.1 什么是CRS

坐标参考系统（Coordinate Reference System，CRS）定义了如何将地球表面的位置映射到二维或三维坐标：

CRS组成
├── 坐标系统
│   ├── 地理坐标系 (经纬度)
│   └── 投影坐标系 (平面坐标)
├── 基准面 (Datum)
│   ├── 参考椭球体
│   └── 大地基准点
└── 投影方法 (仅投影坐标系)
    ├── 投影类型
    └── 投影参数

8.1.2 地理坐标系 vs 投影坐标系

特性 地理坐标系 投影坐标系
坐标单位 度 (°) 米/英尺
形状 球面/椭球面 平面
变形 无 有（面积/角度/距离）
适用范围 全球 特定区域
示例 WGS84 (EPSG:4326) UTM (EPSG:32650)

8.1.3 常用坐标系

全球坐标系：

EPSG代码 名称 说明
4326 WGS 84 GPS坐标系，全球通用
4490 CGCS 2000 中国国家大地坐标系
3857 Web Mercator Web地图标准

中国常用坐标系：

EPSG代码 名称 说明
4490 CGCS2000 2000国家大地坐标系（地理）
4528-4554 CGCS2000 3度带 高斯-克吕格投影
4502-4512 CGCS2000 6度带 高斯-克吕格投影
2432-2442 北京54 3度带 旧坐标系
2401-2414 西安80 3度带 旧坐标系

UTM分带：

带号 EPSG (北半球) 覆盖范围
49N 32649 108°-114° E
50N 32650 114°-120° E
51N 32651 120°-126° E

8.1.4 EPSG代码系统

EPSG (European Petroleum Survey Group) 提供了标准化的CRS编码：

from qgis.core import QgsCoordinateReferenceSystem

# 使用EPSG代码创建CRS
crs = QgsCoordinateReferenceSystem("EPSG:4326")
print(f"名称: {crs.description()}")
print(f"EPSG: {crs.authid()}")
print(f"单位: {crs.mapUnits()}")
print(f"是否地理坐标系: {crs.isGeographic()}")

8.2 QGIS中的CRS管理

8.2.1 项目CRS设置

设置方式：

项目 > 属性 > CRS
或右下角状态栏点击CRS按钮

配置选项：

  • 无CRS（无投影）
  • 预定义CRS
  • 自定义CRS
  • 从图层获取

8.2.2 图层CRS

查看图层CRS：

layer = iface.activeLayer()
crs = layer.crs()
print(f"图层CRS: {crs.authid()}")

设置图层CRS（不转换数据）：

layer.setCrs(QgsCoordinateReferenceSystem("EPSG:4326"))

8.2.3 即时投影（OTF）

QGIS支持即时投影，允许不同CRS的图层同时显示：

启用/禁用：

项目 > 属性 > CRS > 启用"即时"CRS转换

工作原理：

图层1 (EPSG:4326) ─┐
                   ├─> 投影到项目CRS ─> 地图画布显示
图层2 (EPSG:32650) ─┘

8.2.4 CRS选择器

打开方式：

  • 图层属性 > 源 > 设置CRS
  • 数据源管理器加载时
  • 导出数据时

搜索CRS：

  • 按名称搜索
  • 按EPSG代码搜索
  • 从图层获取
  • 从收藏夹选择

8.3 坐标转换

8.3.1 基本坐标转换

from qgis.core import (
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsProject,
    QgsPointXY
)

# 定义源和目标CRS
source_crs = QgsCoordinateReferenceSystem("EPSG:4326")
dest_crs = QgsCoordinateReferenceSystem("EPSG:32650")

# 创建转换对象
transform = QgsCoordinateTransform(
    source_crs,
    dest_crs,
    QgsProject.instance()
)

# 转换点坐标
point_wgs84 = QgsPointXY(116.397, 39.908)  # 北京天安门
point_utm = transform.transform(point_wgs84)
print(f"WGS84: {point_wgs84}")
print(f"UTM50N: {point_utm}")

# 反向转换
reverse_transform = QgsCoordinateTransform(
    dest_crs,
    source_crs,
    QgsProject.instance()
)
point_back = reverse_transform.transform(point_utm)
print(f"转回: {point_back}")

8.3.2 转换范围

from qgis.core import QgsRectangle

# 转换矩形范围
extent_wgs84 = QgsRectangle(116.0, 39.0, 117.0, 40.0)
extent_utm = transform.transformBoundingBox(extent_wgs84)
print(f"原范围: {extent_wgs84}")
print(f"转换后: {extent_utm}")

8.3.3 转换几何对象

from qgis.core import QgsGeometry

# 转换几何
geom_wgs84 = QgsGeometry.fromPointXY(QgsPointXY(116.0, 39.0))
geom_utm = geom_wgs84.transform(transform)
# 注意：transform方法就地修改几何，返回错误码

# 或创建副本转换
geom_copy = QgsGeometry(geom_wgs84)
geom_copy.transform(transform)

8.3.4 批量转换图层

import processing

# 使用Processing重投影
result = processing.run("native:reprojectlayer", {
    'INPUT': input_layer,
    'TARGET_CRS': 'EPSG:32650',
    'OUTPUT': '/path/to/reprojected.shp'
})

8.4 基准面转换

8.4.1 什么是基准面转换

不同大地基准面之间的转换需要特定的转换参数（7参数或3参数）。

常见转换：

  • WGS84 ↔ CGCS2000（差异很小）
  • WGS84 ↔ 北京54
  • WGS84 ↔ 西安80
  • CGCS2000 ↔ 北京54

8.4.2 配置基准面转换

GUI配置：

设置 > 选项 > CRS和变换 > 基准面转换

查看可用转换：

from qgis.core import QgsCoordinateTransform, QgsDatumTransform

source = QgsCoordinateReferenceSystem("EPSG:4326")
dest = QgsCoordinateReferenceSystem("EPSG:4490")

# 获取可用转换
ops = QgsCoordinateTransform.availableTransformations(source, dest)
for op in ops:
    print(f"转换: {op.name}")
    print(f"  精度: {op.accuracy}")
    print(f"  范围: {op.bounds}")

8.4.3 使用特定转换

from qgis.core import QgsCoordinateTransformContext

# 创建转换上下文
context = QgsCoordinateTransformContext()

# 添加特定转换
context.addCoordinateOperation(
    source_crs,
    dest_crs,
    "+proj=pipeline +step +proj=..."  # PROJ字符串
)

# 使用上下文创建转换
transform = QgsCoordinateTransform(source_crs, dest_crs, context)

8.5 自定义坐标系

8.5.1 创建自定义CRS

GUI方式：

设置 > 自定义投影
点击"添加"
输入名称和投影参数

参数格式（PROJ字符串）：

+proj=tmerc +lat_0=0 +lon_0=117 +k=1 +x_0=500000 +y_0=0 +ellps=GRS80 +units=m +no_defs

8.5.2 常用PROJ参数

参数 说明 示例
+proj 投影类型 tmerc, utm, lcc
+lat_0 纬度原点 0
+lon_0 中央经线 117
+k 比例因子 1, 0.9996
+x_0 东偏移 500000
+y_0 北偏移 0
+ellps 椭球体 GRS80, WGS84
+datum 基准面 WGS84
+units 单位 m
+zone UTM带号 50

8.5.3 中国自定义投影示例

CGCS2000 3度带 117°E中央经线：

+proj=tmerc +lat_0=0 +lon_0=117 +k=1 +x_0=39500000 +y_0=0 +ellps=GRS80 +units=m +no_defs

北京54 3度带：

+proj=tmerc +lat_0=0 +lon_0=117 +k=1 +x_0=500000 +y_0=0 +ellps=krass +towgs84=15.8,-154.4,-82.3,0,0,0,0 +units=m +no_defs

8.5.4 PyQGIS创建自定义CRS

from qgis.core import QgsCoordinateReferenceSystem

# 从PROJ字符串创建
proj_string = "+proj=tmerc +lat_0=0 +lon_0=117 +k=1 +x_0=500000 +y_0=0 +ellps=GRS80 +units=m +no_defs"
crs = QgsCoordinateReferenceSystem()
crs.createFromProj(proj_string)

# 从WKT创建
wkt = '''
PROJCS["Custom_TM",
    GEOGCS["GCS_CGCS_2000",
        DATUM["China_2000",
            SPHEROID["CGCS2000",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",117],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
'''
crs = QgsCoordinateReferenceSystem()
crs.createFromWkt(wkt)

# 保存自定义CRS
if crs.isValid():
    crs.saveAsUserCrs("我的自定义投影")

8.6 CRS在数据操作中的应用

8.6.1 加载数据时指定CRS

from qgis.core import QgsVectorLayer

# 方法1：加载后设置CRS（不转换数据）
layer = QgsVectorLayer("/path/to/data.shp", "layer", "ogr")
layer.setCrs(QgsCoordinateReferenceSystem("EPSG:4490"))

# 方法2：通过URI指定
uri = "/path/to/data.shp|crs=EPSG:4490"
layer = QgsVectorLayer(uri, "layer", "ogr")

8.6.2 导出时转换CRS

from qgis.core import QgsVectorFileWriter, QgsCoordinateTransformContext

options = QgsVectorFileWriter.SaveVectorOptions()
options.driverName = "ESRI Shapefile"
options.fileEncoding = "UTF-8"

# 指定目标CRS
dest_crs = QgsCoordinateReferenceSystem("EPSG:32650")
context = QgsCoordinateTransformContext()

error = QgsVectorFileWriter.writeAsVectorFormatV3(
    layer,
    "/path/to/output.shp",
    context,
    options
)

# 使用Processing
result = processing.run("native:reprojectlayer", {
    'INPUT': layer,
    'TARGET_CRS': 'EPSG:32650',
    'OUTPUT': '/path/to/output.shp'
})

8.6.3 空间分析中的CRS

# 确保图层CRS一致
if layer1.crs() != layer2.crs():
    # 重投影
    result = processing.run("native:reprojectlayer", {
        'INPUT': layer2,
        'TARGET_CRS': layer1.crs().authid(),
        'OUTPUT': 'memory:'
    })
    layer2 = result['OUTPUT']

# 然后进行空间分析
result = processing.run("native:intersection", {
    'INPUT': layer1,
    'OVERLAY': layer2,
    'OUTPUT': 'memory:'
})

8.7 坐标计算与测量

8.7.1 距离计算

from qgis.core import QgsDistanceArea

# 创建距离计算器
distance = QgsDistanceArea()
distance.setSourceCrs(layer.crs(), QgsProject.instance().transformContext())
distance.setEllipsoid('GRS80')  # 设置椭球体

# 计算两点距离
point1 = QgsPointXY(116.0, 39.0)
point2 = QgsPointXY(117.0, 40.0)

dist = distance.measureLine(point1, point2)
print(f"距离: {dist} 米")
print(f"距离: {dist/1000} 公里")

# 转换单位
dist_km = distance.convertLengthMeasurement(dist, QgsUnitTypes.DistanceKilometers)

8.7.2 面积计算

# 计算多边形面积
polygon = QgsGeometry.fromWkt("POLYGON((116 39, 117 39, 117 40, 116 40, 116 39))")

area = distance.measureArea(polygon)
print(f"面积: {area} 平方米")

# 转换单位
area_km2 = distance.convertAreaMeasurement(area, QgsUnitTypes.AreaSquareKilometers)
print(f"面积: {area_km2} 平方公里")

8.7.3 周长计算

perimeter = distance.measurePerimeter(polygon)
print(f"周长: {perimeter} 米")

8.7.4 测地线计算

# 计算测地线（大圆弧）
from qgis.core import QgsGeometry, QgsPointXY

# 创建测地线
point1 = QgsPointXY(116.0, 39.0)
point2 = QgsPointXY(-122.0, 37.0)  # 旧金山

# 使用QgsDistanceArea获取测地线距离
dist = distance.measureLine(point1, point2)
print(f"北京到旧金山测地线距离: {dist/1000:.0f} 公里")

8.8 常见问题与解决方案

8.8.1 CRS不匹配

问题：图层显示位置错误
解决：

# 检查图层CRS
print(f"图层CRS: {layer.crs().authid()}")

# 如果CRS错误，设置正确的CRS
layer.setCrs(QgsCoordinateReferenceSystem("EPSG:4326"))

8.8.2 坐标偏移

问题：数据有几十米到几百米的偏移
原因：基准面转换参数不正确
解决：

# 使用正确的基准面转换
context = QgsProject.instance().transformContext()
# 或在设置中配置默认转换

8.8.3 无法识别CRS

问题：显示"未知CRS"
解决：

# 手动设置CRS
layer.setCrs(QgsCoordinateReferenceSystem("EPSG:4326"))

# 或创建.prj文件
# 导出为带CRS信息的格式

8.8.4 跨带问题

问题：数据跨越多个UTM带
解决方案：

  1. 使用地理坐标系存储
  2. 分析时投影到合适的本地坐标系
  3. 使用自定义投影
# 计算最佳UTM带
def get_utm_zone(lon):
    return int((lon + 180) / 6) + 1

def get_utm_epsg(lon, lat):
    zone = get_utm_zone(lon)
    if lat >= 0:
        return f"EPSG:326{zone:02d}"
    else:
        return f"EPSG:327{zone:02d}"

8.8.5 中国坐标偏移（火星坐标）

问题：使用百度/高德地图时坐标偏移
原因：GCJ-02（火星坐标）加密
解决：需要使用专门的转换库

# 提示：这需要第三方库，如 coordTransform_py
# QGIS原生不支持GCJ-02转换

# 示例概念（需安装相关库）
# from coord_convert import wgs84_to_gcj02, gcj02_to_wgs84
# lng_gcj, lat_gcj = wgs84_to_gcj02(lng_wgs, lat_wgs)

8.9 最佳实践

8.9.1 CRS选择建议

场景 推荐CRS 原因
数据存储 CGCS2000 (EPSG:4490) 中国标准
Web显示 Web Mercator (EPSG:3857) Web兼容
距离/面积计算 本地投影坐标系 准确计算
全球数据 WGS84 (EPSG:4326) 通用标准
局部分析 UTM或高斯-克吕格 变形小

8.9.2 工作流程建议

  1. 统一源数据CRS：确保所有输入数据使用相同CRS
  2. 设置正确的项目CRS：选择适合工作区域的投影
  3. 分析前检查CRS：确保所有图层CRS一致
  4. 保存转换参数：记录使用的基准面转换
  5. 导出时注明CRS：确保输出数据包含完整CRS信息

8.9.3 文档记录

# 记录CRS信息
def document_crs_info(layer):
    crs = layer.crs()
    info = {
        "名称": crs.description(),
        "EPSG": crs.authid(),
        "PROJ": crs.toProj(),
        "WKT": crs.toWkt(),
        "单位": str(crs.mapUnits()),
        "是否地理坐标系": crs.isGeographic()
    }
    return info

# 打印信息
for key, value in document_crs_info(layer).items():
    print(f"{key}: {value}")

8.10 小结

本章详细介绍了QGIS中坐标参考系统和投影的使用：

关键要点：

  1. 理解地理坐标系和投影坐标系的区别
  2. 熟悉常用坐标系及其EPSG代码
  3. 掌握坐标转换和基准面转换方法
  4. 能够创建和使用自定义坐标系
  5. 了解CRS在空间分析中的重要性
  6. 掌握常见问题的解决方法

正确使用CRS是所有GIS操作的基础。

上一章：第07章 栅格数据处理

下一章：第09章 地图样式与符号化

posted @ 2026-01-08 14:04  我才是银古  阅读(24)  评论(0)    收藏  举报