How a plot a raster opened through gdal using matplotlib

Question

I have a .tif file that I tried (I did not get an error code... but I am not sure it worked) to georeference with gdal geotransform function using the following code:

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

I then try to plot it:

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

but get the following:

---------------------------------------------------------------------------
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

Could someone help me figure out whether I did the geotransformation properly and how I can plot it?

The original file is an .e00 file called "Bathymetry (depth in meters)" (availablehere) that I had to convert to .tif using arcmap.

Answer 1

Matplotlib expects a Numpy array, not a GDAL Dataset object. So you'll need to read the data from the dataset first (using .ReadAsArray()), before plotting it with 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

Mask the nodata values:

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

Calculate the extent:

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

And plot the result:

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]')