Inaccurate georeferencing when plotting raster files using rasterio

Question

I am trying to develop a plot of a raster image with a background map. I used rasterio for reading the raster file and contextily for obtaining the background map following this tutorial.

The plotted map is wrong in terms of geocoordinates. To elaborate, the raster is not aligning with the feature it is representing (river) webmap. How can I fix this? What is going wrong?

Google colab file with a reproducible example.

pth_raster = '/content/raster_data_resampled.tif'
data, out_path, affine = reproject_raster(in_path=pth_raster, band_to_read=1, out_folder="./raster_temp", destination_crs='epsg:3857') #<-- for webmap epsg:3857
src = rasterio.open(out_path)
src.crs

>CRS.from_epsg(3857)

basemap, basemap_extent = ctx.bounds2img(*src.bounds)

import matplotlib.cm as cm
my_cmap = cm.jet
my_cmap.set_under('k', alpha=0)

plt.figure(figsize=(10,10))
plt.imshow(basemap, extent=basemap_extent)
plt.imshow(src.read(1), extent=basemap_extent, cmap=my_cmap, interpolation='none', clim=[-900, 1e5]) #<-- clim for removing no data (-9999) 
plt.show()

Answer 1

Your basemap_extent is different from src.bounds.

src_extent = src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top

print(basemap_extent)
print(src_extent)

(-10610682.518435026,-10604567.556172213, 3464737.618110468, 3470852.5803732835),
(-10610046.586014193, -10604771.70913406, 3465296.2656639265, 3470571.1425440577)

Try passing src_extent to plt.imshow:

plt.figure(figsize=(10,10))
plt.imshow(basemap, extent=basemap_extent)
plt.imshow(src.read(1), extent=src_extent, cmap=my_cmap, interpolation='none', clim=[-900, 1e5])
plt.show()