0

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

enter image description here

3

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

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.