GDAL hillshade doesn't match with Matplotlib Basemap coastline

Question

I try to plot hillshade image processed with GDAL with matplotib basemap coastline, but the resulting projection did not match. Here are my code:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from osgeo import gdal

bounds = [125.25, 131.5, -4.5, -2.5]
dem = gdal.Open('raster.tif')
hillshade = gdal.DEMProcessing('hillshade.tif', dem, 'hillshade',
                               computeEdges = True,
                              multiDirectional = True, zFactor=2)
hillshadeArray = hillshade.GetRasterBand(1).ReadAsArray()
elevArray = dem.GetRasterBand(1).ReadAsArray()

plt.figure(figsize=(40,15))
m = Basemap(epsg = '4326', 
            resolution = 'h', 
            llcrnrlon = bounds[0], 
            llcrnrlat = bounds[2],
            urcrnrlon = bounds[1], 
            urcrnrlat = bounds[3])
m.imshow(hillshadeArray, cmap='gray', origin='upper')
m.imshow(elev, cmap='terrain', alpha=0.33, origin='upper')
m.drawcoastlines(color='black', linewidth=2)
m.drawmeridians(np.arange(125, 132, 0.25), labels=[0,0,0,1])
m.drawparallels(np.arange(-5, -2, 0.25), labels=[1,0,0,0])
plt.colorbar()
plt.show()

and here are the results:

What is the issue here and how to fix it?

Answer 1

When you read something as array it looses all georeferencing informaton. The position of your maps seems mostly ok, but the stretch in the horizontal direction is not because the cell size in x direction should not be the same as in y direction.

Check raster.tif CRS with gdalinfo raster.tif.

Try using a projected coordinate system e.g. EPSG:3857 for DEM, Hillshade and on basemap.

Another solution might be to read the raster as array, as you did, but with xoff, yoff, xsize, ysize parameters that match the upper left corner of your map extent in pixels, and number of pixels in x and y directons.