第06章 - 坐标系统与投影转换
第06章：坐标系统与投影转换
6.1 坐标系统基础
6.1.1 什么是坐标系统
坐标系统是用于确定地球表面或空间中点位置的参考框架。在地理信息系统中，坐标系统是所有空间数据的基础。
┌─────────────────────────────────────────────────────────────┐
│ 坐标系统分类 │
├─────────────────────────────────────────────────────────────┤
│ │
│ ┌─────────────────────────────────────────────────────┐ │
│ │ 地理坐标系 (Geographic) │ │
│ │ │ │
│ │ • 使用经度和纬度 │ │
│ │ • 单位：度 (degrees) │ │
│ │ • 基于椭球体模型 │ │
│ │ • 示例：WGS84 (EPSG:4326) │ │
│ └─────────────────────────────────────────────────────┘ │
│ │ │
│ │ 投影转换 │
│ ▼ │
│ ┌─────────────────────────────────────────────────────┐ │
│ │ 投影坐标系 (Projected) │ │
│ │ │ │
│ │ • 使用平面坐标 (X, Y) │ │
│ │ • 单位：米 (meters) 或英尺 │ │
│ │ • 基于特定投影方法 │ │
│ │ • 示例：Web Mercator (EPSG:3857) │ │
│ └─────────────────────────────────────────────────────┘ │
│ │
└─────────────────────────────────────────────────────────────┘
6.1.2 地理坐标系组成
地理坐标系由以下要素定义：
| 要素 | 描述 | 示例 |
|---|---|---|
| 椭球体 | 地球形状的数学模型 | WGS84椭球体、GRS80椭球体 |
| 基准面 | 椭球体与地球的定位关系 | WGS84基准面、北京54基准面 |
| 本初子午线 | 经度起算的起始子午线 | 格林威治子午线 |
| 角度单位 | 坐标的测量单位 | 度、弧度 |
常用椭球体参数：
| 椭球体 | 长半轴 (a) | 扁率 (f) |
|---|---|---|
| WGS84 | 6378137.0 | 1/298.257223563 |
| GRS80 | 6378137.0 | 1/298.257222101 |
| 克拉索夫斯基 | 6378245.0 | 1/298.3 |
| CGCS2000 | 6378137.0 | 1/298.257222101 |
6.1.3 投影坐标系组成
投影坐标系在地理坐标系基础上增加了投影参数：
| 要素 | 描述 |
|---|---|
| 基础地理坐标系 | 投影的源坐标系 |
| 投影方法 | 将曲面投影到平面的方法 |
| 投影参数 | 中央经线、标准纬线、原点等 |
| 线性单位 | 米、英尺等 |
| 假东/假北 | 避免负坐标的偏移量 |
6.1.4 常见投影类型
┌─────────────────────────────────────────────────────────────┐
│ 投影类型分类 │
├─────────────────────────────────────────────────────────────┤
│ │
│ 按投影面类型： │
│ ├── 圆柱投影 (Cylindrical) │
│ │ └── 墨卡托、UTM、Web Mercator │
│ ├── 圆锥投影 (Conic) │
│ │ └── 兰伯特等角圆锥、阿尔伯斯等积圆锥 │
│ └── 方位投影 (Azimuthal) │
│ └── 立体投影、正射投影 │
│ │
│ 按保持特性： │
│ ├── 等角投影 (保持形状) │
│ │ └── 墨卡托、兰伯特等角圆锥 │
│ ├── 等积投影 (保持面积) │
│ │ └── 阿尔伯斯等积圆锥、高尔-彼得斯 │
│ └── 等距投影 (保持距离) │
│ └── 等距圆柱投影 │
│ │
└─────────────────────────────────────────────────────────────┘
6.2 EPSG 代码与 WKT
6.2.1 EPSG 代码系统
EPSG（European Petroleum Survey Group）代码是坐标系统的标准标识符：
| EPSG 代码 | 名称 | 类型 | 描述 |
|---|---|---|---|
| 4326 | WGS 84 | 地理 | 全球通用地理坐标系 |
| 4490 | CGCS2000 | 地理 | 中国国家大地坐标系 |
| 3857 | Web Mercator | 投影 | Web地图投影 |
| 32650 | WGS 84 / UTM zone 50N | 投影 | UTM 50带北半球 |
| 4547 | CGCS2000 / 3-degree Gauss-Kruger zone 39 | 投影 | 高斯-克吕格39带 |
| 2154 | RGF93 / Lambert-93 | 投影 | 法国兰伯特投影 |
| 27700 | OSGB 1936 / British National Grid | 投影 | 英国国家网格 |
6.2.2 WKT 格式
WKT（Well-Known Text）是描述坐标系统的文本格式：
WKT1 格式（传统）：
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]]
WKT2 格式（现代，GDAL 3.0+）：
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
6.2.3 PROJ 字符串
PROJ 字符串是另一种简洁的坐标系统定义格式：
# WGS 84
+proj=longlat +datum=WGS84 +no_defs
# Web Mercator
+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
# UTM Zone 50N
+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs
6.3 OSR 模块详解
6.3.1 创建空间参考
from osgeo import osr
# 方法1：从 EPSG 代码创建
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
# 方法2：从 WKT 创建
wkt = """GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]]"""
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
# 方法3：从 PROJ 字符串创建
srs = osr.SpatialReference()
srs.ImportFromProj4('+proj=longlat +datum=WGS84 +no_defs')
# 方法4：从 ESRI WKT 创建
srs = osr.SpatialReference()
srs.ImportFromESRI(['GEOGCS["GCS_WGS_1984"...'])
# 方法5：从 URL 创建
srs = osr.SpatialReference()
srs.ImportFromUrl('http://spatialreference.org/ref/epsg/4326/ogcwkt/')
# 方法6：自定义创建
srs = osr.SpatialReference()
srs.SetGeogCS(
"WGS 84", # 地理坐标系名称
"WGS_1984", # 基准面名称
"WGS 84", # 椭球体名称
6378137, # 长半轴
298.257223563 # 反扁率
)
6.3.2 查询空间参考属性
from osgeo import osr
srs = osr.SpatialReference()
srs.ImportFromEPSG(32650) # UTM Zone 50N
# 基本信息
print(f"名称: {srs.GetName()}")
print(f"是地理坐标系: {srs.IsGeographic()}")
print(f"是投影坐标系: {srs.IsProjected()}")
print(f"是本地坐标系: {srs.IsLocal()}")
# 获取权威信息
print(f"权威名称: {srs.GetAuthorityName(None)}")
print(f"权威代码: {srs.GetAuthorityCode(None)}")
# 获取单位
linear_units = srs.GetLinearUnits()
linear_units_name = srs.GetLinearUnitsName()
print(f"线性单位: {linear_units_name} ({linear_units})")
angular_units = srs.GetAngularUnits()
angular_units_name = srs.GetAngularUnitsName()
print(f"角度单位: {angular_units_name}")
# 获取椭球体参数
semi_major = srs.GetSemiMajor()
semi_minor = srs.GetSemiMinor()
inv_flattening = srs.GetInvFlattening()
print(f"长半轴: {semi_major}")
print(f"短半轴: {semi_minor}")
print(f"反扁率: {inv_flattening}")
# 获取投影参数（如果是投影坐标系）
if srs.IsProjected():
proj_name = srs.GetAttrValue('PROJECTION')
print(f"投影方法: {proj_name}")
# 获取特定参数
central_meridian = srs.GetProjParm(osr.SRS_PP_CENTRAL_MERIDIAN)
false_easting = srs.GetProjParm(osr.SRS_PP_FALSE_EASTING)
false_northing = srs.GetProjParm(osr.SRS_PP_FALSE_NORTHING)
scale_factor = srs.GetProjParm(osr.SRS_PP_SCALE_FACTOR)
print(f"中央经线: {central_meridian}")
print(f"假东: {false_easting}")
print(f"假北: {false_northing}")
print(f"比例因子: {scale_factor}")
6.3.3 导出空间参考
from osgeo import osr
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
# 导出为 WKT
wkt = srs.ExportToWkt()
print("WKT:")
print(wkt)
# 导出为格式化的 WKT
pretty_wkt = srs.ExportToPrettyWkt()
print("\n格式化 WKT:")
print(pretty_wkt)
# 导出为 PROJ 字符串
proj4 = srs.ExportToProj4()
print(f"\nPROJ4: {proj4}")
# 导出为 XML
xml = srs.ExportToXML()
print(f"\nXML (前200字符): {xml[:200]}...")
# 导出为 MICoordSys（MapInfo）
micoordsys = srs.ExportToMICoordSys()
print(f"\nMapInfo: {micoordsys}")
# 导出为 USGS
usgs = srs.ExportToUSGS()
print(f"\nUSGS: {usgs}")
6.3.4 空间参考比较
from osgeo import osr
srs1 = osr.SpatialReference()
srs1.ImportFromEPSG(4326)
srs2 = osr.SpatialReference()
srs2.ImportFromProj4('+proj=longlat +datum=WGS84 +no_defs')
# 严格比较
is_same = srs1.IsSame(srs2)
print(f"严格相同: {is_same}")
# 仅比较地理坐标系部分
is_same_geogcs = srs1.IsSameGeogCS(srs2)
print(f"地理坐标系相同: {is_same_geogcs}")
# 仅比较垂直坐标系部分
is_same_vertcs = srs1.IsSameVertCS(srs2)
print(f"垂直坐标系相同: {is_same_vertcs}")
6.4 坐标转换
6.4.1 基本坐标转换
from osgeo import osr
def transform_point(x, y, src_epsg, dst_epsg):
"""转换单个坐标点"""
# 创建源和目标坐标系
src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(src_epsg)
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(dst_epsg)
# 创建坐标转换对象
transform = osr.CoordinateTransformation(src_srs, dst_srs)
# 执行转换
result = transform.TransformPoint(x, y)
return result[0], result[1]
# WGS84 转 Web Mercator
x, y = 116.4, 39.9
new_x, new_y = transform_point(x, y, 4326, 3857)
print(f"WGS84 ({x}, {y}) -> Web Mercator ({new_x:.2f}, {new_y:.2f})")
# Web Mercator 转 WGS84
x2, y2 = transform_point(new_x, new_y, 3857, 4326)
print(f"Web Mercator ({new_x:.2f}, {new_y:.2f}) -> WGS84 ({x2}, {y2})")
6.4.2 批量坐标转换
from osgeo import osr
import numpy as np
def transform_points(points, src_epsg, dst_epsg):
"""批量转换坐标点"""
src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(src_epsg)
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(dst_epsg)
transform = osr.CoordinateTransformation(src_srs, dst_srs)
# 转换点列表
results = []
for point in points:
if len(point) == 2:
result = transform.TransformPoint(point[0], point[1])
results.append((result[0], result[1]))
elif len(point) == 3:
result = transform.TransformPoint(point[0], point[1], point[2])
results.append((result[0], result[1], result[2]))
return results
def transform_points_array(xs, ys, src_epsg, dst_epsg):
"""使用数组批量转换（更高效）"""
src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(src_epsg)
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(dst_epsg)
transform = osr.CoordinateTransformation(src_srs, dst_srs)
# 使用 TransformPoints 方法（更高效）
points = list(zip(xs, ys))
transformed = transform.TransformPoints(points)
new_xs = [p[0] for p in transformed]
new_ys = [p[1] for p in transformed]
return new_xs, new_ys
# 示例
points = [(116.4, 39.9), (121.5, 31.2), (113.3, 23.1)]
transformed = transform_points(points, 4326, 3857)
for orig, trans in zip(points, transformed):
print(f"{orig} -> ({trans[0]:.2f}, {trans[1]:.2f})")
6.4.3 几何对象转换
from osgeo import ogr, osr
def transform_geometry(geom, src_epsg, dst_epsg):
"""转换几何对象的坐标系"""
src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(src_epsg)
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(dst_epsg)
transform = osr.CoordinateTransformation(src_srs, dst_srs)
# 克隆几何对象（避免修改原对象）
new_geom = geom.Clone()
new_geom.Transform(transform)
return new_geom
# 转换点
point = ogr.CreateGeometryFromWkt('POINT(116.4 39.9)')
point_3857 = transform_geometry(point, 4326, 3857)
print(f"点转换: {point.ExportToWkt()} -> {point_3857.ExportToWkt()}")
# 转换多边形
polygon = ogr.CreateGeometryFromWkt(
'POLYGON((116 39, 117 39, 117 40, 116 40, 116 39))'
)
polygon_3857 = transform_geometry(polygon, 4326, 3857)
print(f"多边形转换完成")
print(f"原面积: {polygon.GetArea()}")
print(f"新面积: {polygon_3857.GetArea()}")
6.4.4 图层整体转换
from osgeo import gdal, ogr, osr
def reproject_layer(src_path, dst_path, dst_epsg, dst_format='ESRI Shapefile'):
"""重投影整个图层"""
options = gdal.VectorTranslateOptions(
format=dst_format,
dstSRS=f'EPSG:{dst_epsg}',
reproject=True
)
ds = gdal.VectorTranslate(dst_path, src_path, options=options)
ds = None
def reproject_layer_manual(src_path, dst_path, dst_epsg):
"""手动重投影图层（更多控制）"""
# 打开源数据
src_ds = ogr.Open(src_path)
src_layer = src_ds.GetLayer()
src_srs = src_layer.GetSpatialRef()
# 创建目标坐标系
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(dst_epsg)
# 创建转换器
transform = osr.CoordinateTransformation(src_srs, dst_srs)
# 创建目标数据源
driver = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = driver.CreateDataSource(dst_path)
dst_layer = dst_ds.CreateLayer(
src_layer.GetName(),
dst_srs,
src_layer.GetGeomType()
)
# 复制字段定义
src_defn = src_layer.GetLayerDefn()
for i in range(src_defn.GetFieldCount()):
dst_layer.CreateField(src_defn.GetFieldDefn(i))
dst_defn = dst_layer.GetLayerDefn()
# 转换每个要素
for src_feature in src_layer:
dst_feature = ogr.Feature(dst_defn)
# 转换几何
geom = src_feature.GetGeometryRef()
if geom:
new_geom = geom.Clone()
new_geom.Transform(transform)
dst_feature.SetGeometry(new_geom)
# 复制属性
for i in range(dst_defn.GetFieldCount()):
dst_feature.SetField(i, src_feature.GetField(i))
dst_layer.CreateFeature(dst_feature)
src_ds = None
dst_ds = None
# 使用示例
reproject_layer('wgs84.shp', 'webmercator.shp', 3857)
6.4.5 栅格重投影
from osgeo import gdal
def reproject_raster(src_path, dst_path, dst_epsg,
resampling='bilinear', resolution=None):
"""重投影栅格数据"""
resample_methods = {
'nearest': gdal.GRA_NearestNeighbour,
'bilinear': gdal.GRA_Bilinear,
'cubic': gdal.GRA_Cubic,
'cubicspline': gdal.GRA_CubicSpline,
'lanczos': gdal.GRA_Lanczos,
'average': gdal.GRA_Average,
'mode': gdal.GRA_Mode,
}
warp_options = {
'dstSRS': f'EPSG:{dst_epsg}',
'resampleAlg': resample_methods.get(resampling, gdal.GRA_Bilinear),
'format': 'GTiff',
'creationOptions': ['COMPRESS=LZW', 'TILED=YES'],
}
if resolution:
warp_options['xRes'] = resolution
warp_options['yRes'] = resolution
ds = gdal.Warp(dst_path, src_path, **warp_options)
ds = None
# 使用示例
reproject_raster('wgs84.tif', 'webmercator.tif', 3857)
reproject_raster('wgs84.tif', 'utm50n.tif', 32650, resolution=30)
6.5 中国常用坐标系统
6.5.1 中国主要坐标系统
| 名称 | EPSG | 类型 | 说明 |
|---|---|---|---|
| 北京54 | 4214 | 地理 | 基于克拉索夫斯基椭球体 |
| 西安80 | 4610 | 地理 | 基于IAG75椭球体 |
| CGCS2000 | 4490 | 地理 | 中国国家大地坐标系 |
| WGS84 | 4326 | 地理 | 全球定位系统坐标系 |
中国投影坐标系：
| 投影类型 | EPSG范围 | 说明 |
|---|---|---|
| 高斯-克吕格3度带 | 4502-4554 | CGCS2000 |
| 高斯-克吕格6度带 | 4491-4501 | CGCS2000 |
| 北京54 3度带 | 2401-2442 | 已弃用 |
| 西安80 3度带 | 2349-2390 | 已弃用 |
6.5.2 坐标系统转换
from osgeo import osr
def create_chinese_srs(name):
"""创建中国常用坐标系"""
srs = osr.SpatialReference()
if name == 'CGCS2000':
srs.ImportFromEPSG(4490)
elif name == 'Beijing54':
srs.ImportFromEPSG(4214)
elif name == 'Xian80':
srs.ImportFromEPSG(4610)
elif name == 'WGS84':
srs.ImportFromEPSG(4326)
elif name.startswith('CGCS2000_GK_'):
# CGCS2000 高斯-克吕格投影
zone = int(name.split('_')[-1])
if zone <= 22:
# 6度带
epsg = 4491 + zone - 13
else:
# 3度带
epsg = 4502 + zone - 25
srs.ImportFromEPSG(epsg)
else:
raise ValueError(f"未知的坐标系: {name}")
return srs
# WGS84 与 CGCS2000 转换
# 注意：WGS84 和 CGCS2000 在厘米级别上非常接近
# 对于大多数应用可以视为相同
wgs84 = create_chinese_srs('WGS84')
cgcs2000 = create_chinese_srs('CGCS2000')
# 检查差异
print(f"WGS84 椭球体: {wgs84.GetAttrValue('SPHEROID')}")
print(f"CGCS2000 椭球体: {cgcs2000.GetAttrValue('SPHEROID')}")
6.5.3 高斯-克吕格投影
from osgeo import osr
def create_gauss_kruger_srs(central_meridian, zone_width=3, datum='CGCS2000'):
"""创建高斯-克吕格投影坐标系"""
srs = osr.SpatialReference()
if datum == 'CGCS2000':
# 设置基础地理坐标系
srs.SetGeogCS(
"China Geodetic Coordinate System 2000",
"China_2000",
"CGCS2000",
6378137.0,
298.257222101
)
elif datum == 'WGS84':
srs.SetGeogCS(
"WGS 84",
"WGS_1984",
"WGS 84",
6378137.0,
298.257223563
)
# 计算带号和假东
if zone_width == 3:
zone = int(central_meridian / 3)
false_easting = zone * 1000000 + 500000
else: # 6度带
zone = int(central_meridian / 6) + 1
false_easting = zone * 1000000 + 500000
# 设置横轴墨卡托投影参数
srs.SetTM(
0, # 原点纬度
central_meridian, # 中央经线
1.0, # 比例因子
false_easting, # 假东
0 # 假北
)
srs.SetLinearUnits("metre", 1.0)
return srs
# 创建适用于北京的 3 度带投影 (中央经线 117°)
beijing_gk = create_gauss_kruger_srs(117, zone_width=3)
print(f"北京高斯-克吕格投影 WKT:")
print(beijing_gk.ExportToPrettyWkt())
6.5.4 Web 墨卡托与 WGS84 转换
from osgeo import osr
import math
def wgs84_to_web_mercator(lon, lat):
"""WGS84 转 Web Mercator (手动计算)"""
x = lon * 20037508.34 / 180
y = math.log(math.tan((90 + lat) * math.pi / 360)) / (math.pi / 180)
y = y * 20037508.34 / 180
return x, y
def web_mercator_to_wgs84(x, y):
"""Web Mercator 转 WGS84 (手动计算)"""
lon = x / 20037508.34 * 180
lat = y / 20037508.34 * 180
lat = 180 / math.pi * (2 * math.atan(math.exp(lat * math.pi / 180)) - math.pi / 2)
return lon, lat
# 使用 GDAL/OSR 转换（更准确）
def convert_wgs84_webmercator(lon, lat, to_webmercator=True):
"""使用 OSR 进行 WGS84 和 Web Mercator 互转"""
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
web_mercator = osr.SpatialReference()
web_mercator.ImportFromEPSG(3857)
if to_webmercator:
transform = osr.CoordinateTransformation(wgs84, web_mercator)
else:
transform = osr.CoordinateTransformation(web_mercator, wgs84)
result = transform.TransformPoint(lon, lat)
return result[0], result[1]
# 测试
lon, lat = 116.4, 39.9
x, y = convert_wgs84_webmercator(lon, lat, to_webmercator=True)
print(f"WGS84 ({lon}, {lat}) -> Web Mercator ({x:.2f}, {y:.2f})")
lon2, lat2 = convert_wgs84_webmercator(x, y, to_webmercator=False)
print(f"Web Mercator ({x:.2f}, {y:.2f}) -> WGS84 ({lon2:.6f}, {lat2:.6f})")
6.6 坐标转换的精度问题
6.6.1 坐标转换误差来源
| 误差来源 | 描述 | 典型量级 |
|---|---|---|
| 基准面差异 | 不同椭球体和定向 | 米到数十米 |
| 投影变形 | 地球曲面到平面的变形 | 取决于位置和投影 |
| 数值精度 | 计算过程中的舍入误差 | 毫米级 |
| 参数精度 | 转换参数的精度 | 厘米到分米 |
6.6.2 七参数转换
from osgeo import osr
def create_seven_param_transform(src_srs, dst_srs,
dx, dy, dz, rx, ry, rz, ds):
"""
创建七参数转换
参数:
dx, dy, dz: 平移参数 (米)
rx, ry, rz: 旋转参数 (弧秒)
ds: 尺度因子 (ppm)
"""
# GDAL 3.0+ 支持通过 TOWGS84 或 CreateCoordinateTransformation 设置
# 这里展示如何在 WKT 中设置 TOWGS84 参数
wkt = f"""
GEOGCS["Custom",
DATUM["Custom_Datum",
SPHEROID["GRS 1980",6378137,298.257222101],
TOWGS84[{dx},{dy},{dz},{rx},{ry},{rz},{ds}]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]]
"""
custom_srs = osr.SpatialReference()
custom_srs.ImportFromWkt(wkt)
return custom_srs
# 北京54 转 WGS84 的参考七参数（示例值，实际应使用测量值）
# beijing54_to_wgs84 = create_seven_param_transform(
# None, None,
# dx=-24.0, dy=123.0, dz=94.0, # 平移
# rx=-0.02, ry=0.25, rz=0.13, # 旋转
# ds=1.1 # 尺度
# )
6.6.3 使用网格文件转换
from osgeo import osr, gdal
def transform_with_grid(src_srs, dst_srs, grid_file):
"""使用网格文件进行高精度转换"""
# 设置 PROJ 网格文件路径
gdal.SetConfigOption('PROJ_LIB', '/path/to/proj/data')
# 某些高精度转换需要下载网格文件
# 例如 NAD27 到 NAD83 的 NADCON 网格
transform = osr.CoordinateTransformation(src_srs, dst_srs)
return transform
# PROJ 网格文件通常存放在 PROJ_LIB 目录下
# 中国地区的高精度转换可能需要专门的网格文件
6.7 实际应用案例
6.7.1 GPS 数据转换
from osgeo import ogr, osr
def process_gps_track(input_file, output_file, output_epsg=32650):
"""
处理 GPS 轨迹数据
- 输入：WGS84 坐标的 GPX/GeoJSON
- 输出：投影坐标系下的 Shapefile
"""
# 打开输入文件
src_ds = ogr.Open(input_file)
src_layer = src_ds.GetLayer()
# 创建坐标转换
src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(4326)
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(output_epsg)
transform = osr.CoordinateTransformation(src_srs, dst_srs)
# 创建输出
driver = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = driver.CreateDataSource(output_file)
dst_layer = dst_ds.CreateLayer('track', dst_srs, ogr.wkbLineString)
# 添加字段
dst_layer.CreateField(ogr.FieldDefn('length_m', ogr.OFTReal))
# 处理每个要素
for src_feature in src_layer:
geom = src_feature.GetGeometryRef()
if geom:
# 转换坐标
new_geom = geom.Clone()
new_geom.Transform(transform)
# 创建要素
dst_feature = ogr.Feature(dst_layer.GetLayerDefn())
dst_feature.SetGeometry(new_geom)
dst_feature.SetField('length_m', new_geom.Length())
dst_layer.CreateFeature(dst_feature)
src_ds = None
dst_ds = None
print(f"GPS 轨迹已转换并保存到 {output_file}")
# 使用示例
# process_gps_track('track.gpx', 'track_utm.shp', 32650)
6.7.2 Web 地图瓦片坐标转换
import math
def tile_to_bounds(z, x, y, srs='WGS84'):
"""
将瓦片坐标转换为边界范围
参数:
z: 缩放级别
x, y: 瓦片坐标
srs: 输出坐标系 ('WGS84' 或 'WebMercator')
"""
n = 2.0 ** z
# 计算 WGS84 边界
lon_min = x / n * 360.0 - 180.0
lon_max = (x + 1) / n * 360.0 - 180.0
lat_max_rad = math.atan(math.sinh(math.pi * (1 - 2 * y / n)))
lat_max = math.degrees(lat_max_rad)
lat_min_rad = math.atan(math.sinh(math.pi * (1 - 2 * (y + 1) / n)))
lat_min = math.degrees(lat_min_rad)
if srs == 'WGS84':
return (lon_min, lat_min, lon_max, lat_max)
elif srs == 'WebMercator':
from osgeo import osr
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
web_mercator = osr.SpatialReference()
web_mercator.ImportFromEPSG(3857)
transform = osr.CoordinateTransformation(wgs84, web_mercator)
min_point = transform.TransformPoint(lon_min, lat_min)
max_point = transform.TransformPoint(lon_max, lat_max)
return (min_point[0], min_point[1], max_point[0], max_point[1])
def wgs84_to_tile(lon, lat, z):
"""将 WGS84 坐标转换为瓦片坐标"""
n = 2.0 ** z
x = int((lon + 180.0) / 360.0 * n)
lat_rad = math.radians(lat)
y = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
return x, y
# 示例
z, x, y = 14, 13533, 6173
bounds = tile_to_bounds(z, x, y)
print(f"瓦片 ({z}/{x}/{y}) 的 WGS84 边界: {bounds}")
lon, lat = 116.4, 39.9
tile_x, tile_y = wgs84_to_tile(lon, lat, 14)
print(f"坐标 ({lon}, {lat}) 在 z=14 时的瓦片: ({tile_x}, {tile_y})")
6.8 本章小结
本章详细介绍了坐标系统与投影转换：
-
坐标系统基础：
- 地理坐标系和投影坐标系的区别
- 椭球体、基准面、投影方法等概念
- 常见投影类型及其特点
-
坐标系统标识：
- EPSG 代码系统
- WKT 格式（WKT1 和 WKT2）
- PROJ 字符串
-
OSR 模块：
- 创建空间参考对象
- 查询和导出空间参考
- 空间参考比较
-
坐标转换：
- 点坐标转换
- 几何对象转换
- 图层和栅格转换
-
中国常用坐标系统：
- CGCS2000、北京54、西安80
- 高斯-克吕格投影
-
精度问题：
- 误差来源
- 七参数转换
- 网格文件转换
6.9 思考与练习
- 解释地理坐标系和投影坐标系的主要区别。
- 为什么 Web Mercator 不适合用于面积计算？
- WGS84 和 CGCS2000 在什么情况下可以视为相同？
- 编写代码计算某点在不同投影下的坐标变化。
- 如何选择适合中国某个省份的投影坐标系？
- 实现一个函数，将整个 Shapefile 文件从一个坐标系转换到另一个。
- 解释七参数转换的原理和适用场景。
浙公网安备 33010602011771号