使用Geopandas、Rasterio和Contextily时栅格与Shapefile无法对齐

我正在完成实验作业，需用Python实现DEM栅格与Shapefile的对齐，后续还要从栅格和多边形图层提取数据至点图层。我熟悉ArcGIS手动操作，但作业要求使用Python或R（我选择Python）。

课程资料显示两者均采用EPSG 3847坐标系，但Shapefile缺失CRS，已通过Geopandas为其设置EPSG 3847；DEM实际为EPSG 3006，我尝试将DEM转换为EPSG 3847，或把Shapefile转换为EPSG 3006，均无法使两者对齐。单独绘制DEM可正常显示，但与Shapefile同图时DEM不出现，而在ArcGIS中两者可正常对齐。

导入依赖库 #

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  

数据读取与坐标系处理 #

# 读取数据文件
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()  

多面板绘图代码 #

### 创建绘图网格  
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)  

### 图1 - 仅底图  
abisveg_3006.plot(ax = ax1, color = 'none')  
cx.add_basemap(ax1, crs = 3006)  
ax1.set_aspect('equal')  
ax1.set_title("AOI底图")  


### 图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('AOI数字高程模型')  


### 图3 - 植被类型  
abisveg_3006.plot(ax = ax3, column = "VEGKOD", cmap = "viridis")  
cx.add_basemap(ax3, crs = 3006)  
ax3.set_aspect('equal')  
ax3.set_title("植被类型")  

单独绘制DEM的代码 #

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

相关截图 #

• 3面板地图缺失DEM：
• 单独DEM显示：

内容的提问来源于stack exchange，提问作者MrBlueSky