如何使用 matplotlib 通过 gdal 打开栅格图

Question

我有一个 .tif 文件，我尝试（我没有收到错误代码......但我不确定它是否有效）使用以下代码使用 gdal geotransform 函数进行地理配准：

raster = gdal.Open("drive/My Drive/raster.tif")
geotransform = raster.GetGeoTransform()
gt = list(geotransform)
gt[0] = -71.9002
gt[3]  = 41.8738
gt[2] = 0
gt[4] = 0
gt[1] = 50
gt[5] =-50

然后我尝试绘制它：

fig, ax = plt.subplots(figsize = (10,10))
rasterio.plot.show(raster, ax=ax)
plt.show()

但得到以下内容：

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-70-d9daf3e89a18> in <module>()
      1 minx, miny, maxx, maxy = ri.geometry.total_bounds
      2 fig, ax = plt.subplots(figsize = (10,10))
----> 3 rasterio.plot.show(raster, ax=ax)
      4 plt.show()

5 frames
/usr/local/lib/python3.7/dist-packages/matplotlib/image.py in set_data(self, A)
    692                 not np.can_cast(self._A.dtype, float, "same_kind")):
    693             raise TypeError("Image data of dtype {} cannot be converted to "
--> 694                             "float".format(self._A.dtype))
    695 
    696         if not (self._A.ndim == 2

TypeError: Image data of dtype object cannot be converted to float

谁能帮我弄清楚我是否正确地进行了地理变换以及如何绘制它？

原始文件是一个名为“Bathymetry（以米为单位的深度）”的 .e00 文件（可用here），我必须使用 arcmap 将其转换为 .tif。

Answer 1

Matplotlib 需要一个 Numpy 数组，而不是一个 GDAL 数据集对象。因此，您需要先从数据集中读取数据（使用 .ReadAsArray()），然后再使用 Matplotlib 对其进行绘图。

infile = r'/vsizip/C:\Temp\bathygridm.zip/bathygridm.e00'

ds = gdal.OpenEx(infile)
gt = ds.GetGeoTransform()
nodata = ds.GetRasterBand(1).GetNoDataValue()
data = ds.ReadAsArray()
ds = None

屏蔽无数据值：

data = np.ma.masked_values(data, nodata)

计算范围：

ys, xs = data.shape
ulx, xres, _, uly, _, yres = gt
extent = [ulx, ulx+xres*xs, uly, uly+yres*ys]

并绘制结果：

fig, ax = plt.subplots(figsize=(5,6), constrained_layout=True, facecolor='w', dpi=86)

cmap = mpl.cm.get_cmap("viridis").copy() # cmap = plt.cm.viridis
cmap.set_bad('#dddddd')

im = ax.imshow(data, extent=extent, cmap=cmap)
cb = fig.colorbar(im, shrink=.5)
cb.set_label('Bathymetry [m]')

如何使用 matplotlib 通过 gdal 打开栅格图

How a plot a raster opened through gdal using matplotlib

python

raster

matplotlib

gdal

geopandas