Raster and Shapefiles not lining up using Geopandas, Rasterio, and Contextily

Question

I am trying to get a DEM raster to line up with a shapefile in Python, but it will not show up no matter what I do. This is for lab exercise, the entire rest of the exercise relies on these lining up, as I will be extracting data from the raster and polygon layers to a point layer.

I know how to do all this "by hand" in ArcGIS, but the point of the exercise is to use R or Python (the professor did an example with R, but we can use whichever, and I have been learning Python the past couple of months for a work project). In the class notes, he says that both files are in EPSG 3847, but the shapefile was missing the CRS, so I added the CRS to it in geopandas.

The DEM appears to be EPSG 3006 (even though it was supposed to be in 3847), so I tried converting it to EPSG 3847 and it still does not show up. So then I tried going the other way and converting the shapefile to EPSG 3006, which did not help either.

import contextily as cx  
import geopandas as gpd  
import rasterio  
from rasterio.plot import show  
from rasterio.crs import CRS  
from rasterio.plot import show as rioshow  
import matplotlib.pyplot as plt  

#data files
abisveg = gpd.read_file(r'/content/drive/MyDrive/Stackoverflow/Sweden/abisveg_polygon.shp')  
abisveg_3847 = abisveg.set_crs(epsg = 3847)  

abisveg_3006 = abisveg_3847.to_crs(epsg = 3006)  

src = rasterio.open(r'/content/drive/MyDrive/Stackoverflow/Sweden/nh_75_6.tif')  
DEM = src.read()  

### creating plot grid  
fig = plt.figure(figsize = (20,20), constrained_layout = True)  
gs = fig.add_gridspec(1,3)  
ax1 = fig.add_subplot(gs[0,0])  
ax2 = fig.add_subplot(gs[0,1], sharex = ax1, sharey = ax1)  
ax3 = fig.add_subplot(gs[0,2], sharex = ax1, sharey = ax1)  

### Plot 1 - Basemap Only  
abisveg_3006.plot(ax = ax1, color = 'none')  
cx.add_basemap(ax1, crs = 3006)  
ax1.set_aspect('equal')  
ax1.set_title("Basemap of AOI")  


### Plot 2 - DEM  
# abisveg_3847.plot(ax = ax2, color = 'none')  
show(DEM, ax=ax2, cmap = "Greys")  
cx.add_basemap(ax2, crs = 3006)  
ax2.set_aspect('equal')  
ax2.set_title('Digitial Elevation Model of AOI')  


### Plot 3 - Vegetation Types  
abisveg_3006.plot(ax = ax3, column = "VEGKOD", cmap = "viridis")  
cx.add_basemap(ax3, crs = 3006)  
ax3.set_aspect('equal')  
ax3.set_title("Vegetation Types")  

3 Panel map with missing DEM: https://i.stack.imgur.com/aIKr8.jpg

Trying to plot the files in Matplotlib has not worked, b/c they do not align at all. I am using contextily for the basemap, and have set the basemap CRS to EPSG 3847 (or 3006, depending on which version of the GIS files I was using). The shapefile shows up in the correct location no matter the projection, but the Raster does not show up. What's weird is that if I open everything up in ArcGIS, it all lines up correctly.

If I plot just the DEM all by itself, it shows up, though I don't know where on the earth it is plotting.

fig = plt.figure(figsize = (10,10), constrained_layout = True)  
show(DEM, cmap = "Greys")  

DEM just by itself: https://i.stack.imgur.com/nQOXj.jpg

I have my code in a colab notebook here:

https://colab.research.google.com/drive/1VAZ3dgf0QS2PPBOl8KJ2FXtB2oRj0qJ8?usp=share_link

The files are here:

https://drive.google.com/drive/folders/1t-xvpIcLOIR9uYXOguJ7KyKqt7wuYSNc?usp=share_link

Answer 1

You could give EOmaps a try... it uses matplotlib/cartopy for plotting and handles re-projecting the data and shapes to the plot-crs

from pathlib import Path
from eomaps import Maps
import geopandas as gpd

p = Path(r"path to the data folder")
# read shapefile
abisveg = gpd.read_file(p / 'abisveg_polygon.shp').set_crs(epsg = 3847)

# create a map in epsg=3006
m = Maps(crs=3006, figsize=(10, 8))
# add stamen-terrain basemap
m.add_wms.OpenStreetMap.add_layer.stamen_terrain()
# plot shapefile (zorder=2 to be on top of the DEM)
m.add_gdf(abisveg, column=abisveg.VEGKOD, cmap="viridis", ec="k", lw=0.2, alpha=0.5, zorder=2)
# plot DEM
m2 = m.new_layer_from_file.GeoTIFF(p / "nh_75_6.tif", cmap="Greys", zorder=1)

m.ax.set_extent((589913.0408156103, 713614.6619114348, 7495264.310799116, 7618965.93189494),
                Maps.CRS.epsg(3006))