光栅和 Shapefile 没有使用 Geopandas、Rasterio 和 Contextily 对齐答案

【问题标题】：Raster and Shapefiles not lining up using Geopandas, Rasterio, and Contextily光栅和 Shapefile 没有使用 Geopandas、Rasterio 和 Contextily 对齐
【发布时间】：2022-11-14 18:44:42
【问题描述】：

我正在尝试让 DEM 栅格与 Python 中的 shapefile 对齐，但无论我做什么，它都不会显示。这是用于实验室练习的，整个练习的其余部分都依赖于这些排列，因为我将从栅格和多边形图层中提取数据到点图层。

我知道如何在 ArcGIS 中“手动”完成所有这些操作，但练习的重点是使用 R 或 Python（教授用 R 做了一个示例，但我们可以使用任何一个，而且我在过去的一对中一直在学习 Python几个月的工作项目）。在课堂笔记中，他说这两个文件都在 EPSG 3847 中，但是 shapefile 缺少 CRS，所以我在 geopandas 中添加了 CRS。

DEM 似乎是 EPSG 3006（即使它应该在 3847 中），所以我尝试将其转换为 EPSG 3847，但它仍然没有出现。所以然后我尝试另一种方式并将shapefile转换为EPSG 3006，这也没有帮助。

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 缺少 DEM 的面板图： https://i.imgur.com/taG2U9Q.jpg

试图在 Matplotlib 中绘制文件没有奏效，因为它们根本不对齐。我在上下文中使用底图，并将底图 CRS 设置为 EPSG 3847（或 3006，取决于我使用的 GIS 文件的版本）。无论投影如何，shapefile 都会显示在正确的位置，但 Raster 不会显示。奇怪的是，如果我在 ArcGIS 中打开所有内容，一切都会正确排列。

如果我只单独绘制 DEM，它就会显示出来，尽管我不知道它在地球上的什么地方绘制。

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

DEM 本身： https://i.imgur.com/KyYu7jc.jpg

我在 colab 笔记本中有我的代码：

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

文件在这里：

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

【问题讨论】：

在show(DEM, ax=ax2, cmap = "Greys") 中，您需要指定适当的选项以将栅格放置在正确的位置、比例和 CRS。

标签： python matplotlib geopandas rasterio epsg

【解决方案1】：

你可以试试 EOmaps...它使用 matplotlib/cartopy 进行绘图并处理将数据和形状重新投影到 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))

【讨论】：