通过GDAL打开的栅格如何用matplotlib绘制？地理变换及绘图报错如何解决？

一、地理变换操作存在的问题 #

你当前的地理配准操作未生效，核心问题是修改完地理变换参数后没有写入GDAL数据集对象：

• 你修改了gt列表的参数，但没有执行raster.SetGeoTransform(tuple(gt))，修改只停留在Python列表层面，没有同步到GDAL打开的栅格对象中，配准不会生效
• 额外注意参数匹配问题：你设置的像素分辨率为50/-50，需和栅格的坐标系单位匹配：如果是地理坐标系（单位为度），50的分辨率完全不符合实际业务场景；如果是投影坐标系（单位为米），需确认你设置的左上角坐标-71.9002/41.8738是对应投影坐标系下的坐标，而非WGS84经纬度

二、报错原因与解决方案 #

报错根因 #

你混用了GDAL和rasterio的数据集对象：rasterio.plot.show仅支持rasterio自身打开的栅格对象、numpy数组两类输入，传入GDAL的数据集对象时无法正确解析栅格数组，返回object类型数据，触发类型转换报错。

解决步骤 #

两种方案二选一即可：

方案1：全程用GDAL读取数据+matplotlib绘图

修改代码如下：

from osgeo import gdal
import matplotlib.pyplot as plt

# 打开栅格并修改地理变换
raster = gdal.Open("drive/My Drive/raster.tif")
geotransform = list(raster.GetGeoTransform())
geotransform[0] = -71.9002
geotransform[3] = 41.8738
geotransform[2] = 0
geotransform[4] = 0
geotransform[1] = 50
geotransform[5] = -50
# 写入修改后的地理变换，这一步是配准生效的核心
raster.SetGeoTransform(tuple(geotransform))

# 读取栅格数组
arr = raster.ReadAsArray()
# 计算地理范围用于显示坐标
x_min = geotransform[0]
x_max = geotransform[0] + geotransform[1] * arr.shape[1]
y_max = geotransform[3]
y_min = geotransform[3] + geotransform[5] * arr.shape[0]

# 绘图
fig, ax = plt.subplots(figsize=(10,10))
ax.imshow(arr, extent=(x_min, x_max, y_min, y_max))
plt.show()

# 如果需要永久保存配准结果，执行以下代码将修改写入文件
raster = None

方案2：全程用rasterio处理

修改代码如下：

import rasterio
from rasterio.plot import show
import matplotlib.pyplot as plt
from affine import Affine

with rasterio.open("drive/My Drive/raster.tif", "r+") as src:
    # 修改地理变换，rasterio用Affine对象存储变换参数
    new_transform = Affine(50, 0, -71.9002, 0, -50, 41.8738)
    src.transform = new_transform
    # 绘图
    fig, ax = plt.subplots(figsize=(10,10))
    show(src, ax=ax)
    plt.show()

内容的提问来源于stack exchange，提问作者nalfahel