光栅和形状文件不使用Geopandas、Rasterio和Contextily排列

Question

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

我知道如何在ArcGIS中“手工”完成所有这些操作，但练习的重点是使用R或Python (教授用R做了一个示例，但我们可以使用哪种方法，过去几个月我一直在为一个工作项目学习Python )。在课堂笔记中，他说这两个文件都在EPSG 3847中，但是shapefile缺少CRS，所以我在地质公园中添加了CRS。

DEM似乎是EPSG 3006 (尽管它应该在3847)，所以我尝试将它转换为EPSG 3847，但它仍然没有出现。于是，我试着用另一种方式将shapefile转换为EPSG3006，这也没有帮助。

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中的文件没有工作，b/c它们根本不对齐。我正在为basemap使用上下文，并将basemap 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

Answer 1

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

​

​