I've been stuck for a few days trying to plot a TIF file over a basemap (Google Maps or Mapbox) using Python. I don't know if this is actually possible because I didn't find examples and nothing very clear about it, the most I found is how to plot vector files (shp) on top of basemaps, but nothing about raster.
I have a raster (rainfall_clipped.tif) that is in EPSG:4326 and I'm trying to plot it on top of a MapBox map using the Contextily package. At first I suspected that it could be a projection problem, but even doing the conversions to EPSG:3857, it didn't work. In the documentation Contextily says that it accepts both projections.
Contextily seems to understand the projections, even because it generates the map with the correct axes, the big problem is that it doesn't plot the TIF file on the map. And it's also reading the raster because it loaded the colorbar with the correct values. I don't know if the layer order is wrong or what else it could be.
Below is the code I am using and then the image that is being generated:
data = rasterio.open("rainfall_clipped.tif")
# Read the bounds
left, bottom, right, top = data.bounds
# Create a figure and axes
fig, ax = plt.subplots(figsize=(10, 10))
# Add the raster data to the plot using imshow
im = ax.imshow(data.read(1), extent=[left, right, bottom, top], cmap="Blues")
# Add a colorbar
fig.colorbar(im, ax=ax)
# Add a basemap using contextily
ctx.add_basemap(ax, crs=data.crs, source=ctx.providers.MapBox(accessToken="my_key", id="mapbox/satellite-v9"))
# Show the plot
plt.show()
This code is generating this image: Output image
What I'm trying to do is something like below: Expected image
Any suggestions what I can do?