Как открыть растр через gdal с помощью matplotlib [закрыто]

avatar
nalfahel
8 августа 2021 в 22:06
424
1
-1

У меня есть файл .tif, который я попытался (я не получил код ошибки... но я не уверен, что это сработало) выполнить географическую привязку с помощью функции геотрансформации gdal, используя следующий код:

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

Может ли кто-нибудь помочь мне понять, правильно ли я выполнил геотрансформацию и как ее построить?

Исходный файл представляет собой файл .e00 под названием «Батиметрия (глубина в метрах)» (доступенздесь), который мне пришлось преобразовать в .tif с помощью arcmap.

Источник

Ответы (1)

avatar
Rutger Kassies
9 августа 2021 в 09:33
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]')

enter image description here

nalfahel
9 августа 2021 в 12:35
0

Привет! Большое спасибо за публикацию ответа! У меня только один вопрос.... диапазон значения растра должен быть от -2 до 65 (или что-то в этом роде) почему цветная полоса начинается с 0? Я смог построить другой растр (те же данные, но большей протяженности из другого источника), и при импорте и чтении его с использованием растра масштаб менялся от -3 до 216 в фактических данных до 0-150 на графике...

Rutger Kassies
10 августа 2021 в 06:31
0

@nalfahel, ты прав, у него должно быть минимум -2. Не уверен, почему мое изображение этого не показывало, сначала я просто обрезал его, добавив vmin=0 в вызов imshow(), прежде чем правильно замаскировать значение nodata. Возможно, я не сохранил свою последнюю версию. Я обновил изображение в посте выше.