Python+GDAL 栅格坐标系转换（自动计算输出像元大小）

news2024/11/18 9:01:36

GDAL对栅格进行坐标系转换不难，直接用gdal.Warp()就可以了

gdal.Warp("output", "input", dstSRS='EPSG:***', xRes=**, yRes=**, targetAlignedPixels=True)

麻烦的是，需要的参数xRes和yRes，gdal.Warp()不能自动计算。坐标系转换后的最佳分辨率的计算方法，在ArcGIS中有描述，下滑最后有原文。此处附链接：https://desktop.arcgis.com/zh-cn/arcmap/latest/tools/environments/how-the-cell-size-projection-method-environment-setting-works.htm
计算距离考虑两种情况：地理坐标系的距离和空间直角坐标系的距离。
地理坐标系的距离计算应计算测地线距离（椭球面距离），空间直角坐标系可直接计算欧氏距离，距离计算的脚本如下：

import numpy as np
from geopy import distance
# 计算测地线距离（输入两组经纬度坐标）
a = distance.distance((110.5, 39.2), (110.4, 38.9))
# 计算欧式距离（输入两组投影坐标）
b = np.sqrt(np.sum((np.array([594892.5,3947407.5]) - np.array([822307.5,3715792.5])) ** 2))

根据ArcGIS的转换单位方法计算像元大小，需要计算原始范围和投影范围的六边距离（六边=四条边+两条对角线），下面是一个通过四至计算六边距离的方法：

from geopy import distance
from typing import List
# 计算六边距离
def calculate_distance(bounds:List[float]):
    [orX,orY,edX,edY] = bounds
    if -90.<= orY <= 90. and -90.<= edY <= 90.:
        a = distance.distance((orX, orY), (edX, edY))
        b = distance.distance((orX, edY), (edX, orY))
        c = distance.distance((orX, orY), (orX, edY))
        d = distance.distance((orX, orY), (edX, orY))
        e = distance.distance((edX, edY), (edX, orY))
        f = distance.distance((edX, edY), (orX, edY))
    else:
        a = np.sqrt(np.sum((np.array([orX, orY]) - np.array([edX, edY])) ** 2))
        b = np.sqrt(np.sum((np.array([orX, edY]) - np.array([edX, orY])) ** 2))
        c = np.sqrt(np.sum((np.array([orX, orY]) - np.array([orX, edY])) ** 2))
        d = np.sqrt(np.sum((np.array([orX, orY]) - np.array([edX, orY])) ** 2))
        e = np.sqrt(np.sum((np.array([edX, edY]) - np.array([edX, orY])) ** 2))
        f = np.sqrt(np.sum((np.array([edX, edY]) - np.array([orX, edY])) ** 2))

    return [a,b,c,d,e,f]

有了上述的计算六边距离，就可以计算输出栅格的像元大小xRes和yRes了，完整的栅格坐标系转换方法如下：

from osgeo import osr,gdal
# 打开栅格并读取属性
raster = gdal.Open(in_file)
cols = raster.RasterXSize
rows = raster.RasterYSize
orX, psX, _, orY, _, psY = raster.GetGeoTransform()
# 计算右上角坐标
edX = orX + psX * cols
edY = orY + psY * rows
# print(orX,orY,edX,edY)
# 计算四边与对角线边长
dist0 = calculate_distance([orX,orY,edX,edY])
# 定义投影转换（转为WGS84）
proj = osr.SpatialReference()
proj.ImportFromWkt(raster.GetProjection())
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
tx = osr.CoordinateTransformation(proj, wgs84)
# 计算投影后四至
(ulx, uly, _) = tx.TransformPoint(orX, orY)
(lrx, lry, _) = tx.TransformPoint(edX, edY)
# 计算四边与对角线边长
dist1 = calculate_distance([ulx,uly,lrx,lry])
# 计算scale
scale = (np.array(dist1)/np.array(dist0)).sum() / 6.
# print(scale)
# 计算投影后分辨率
reX = psX * scale
reY = -psY * scale
# 重投影
gdal.Warp(out_file, in_file, dstSRS='EPSG:4326', xRes=reX, yRes=reY, targetAlignedPixels=True)