3

I have a DEM and a polygon, both in the NAD83 Datum. I need to plot both at the same time.

Plotting the polygon

import os
import rasterio
import geopandas as gpd
import contextily as cx
import matplotlib.pyplot as plt
from osgeo import ogr2ogr

site_bdry = gpd.read_file(boundary_json)
site_bdry = site_bdry.to_crs("EPSG:26913") # NAD83 / UTM zone 13N
ax = site_bdry.plot(figsize=(10,11), alpha=0)
site_bdry.plot(ax=ax, facecolor="none", edgecolor="red", linewidth=2)
cx.add_basemap(ax, crs=site_bdry.crs.to_string())

enter image description here

Plotting the DEM

dem = rasterio.open(dem_path)
ax = plt.subplots(figsize=(20, 10))
plt.imshow(dem.read(1), cmap="Greys")
plt.show()

enter image description here

Plotting the two together like this Plot shapefile on top of raster using plot and imshow from matplotlib and this How can I superimpose a geopandas dataframe on a raster plot?

site_bdry = gpd.read_file(boundary_json)
site_bdry = site_bdry.to_crs("EPSG:26913")
site_bdry_buff = site_bdry["geometry"]
dem = rasterio.open(dem_path)
fig, ax = plt.subplots(figsize=(5, 6))
plt.imshow(dem.read(1), cmap="Greys")
site_bdry.plot(ax=ax, facecolor="none", edgecolor="red", linewidth=2)
plt.show()

enter image description here

0

2 Answers 2

0

Compare the two plotted figures. You can see that the numbers on the x-axis and y-axis are completely different. For raster the proper crs was not entered. If your DEM's coordinates are EPSG:26913

from rasterio.crs import CRS    

dem = rasterio.open(dem_path, 'r+')
dem.crs = CRS.from_epsg(26913)

If the crs of the DEM are different, you have to go through the warp process. Find related examples.

0

Try rasterio.plot show:

Rasterio also provides rasterio.plot.show() to perform common tasks such as displaying multi-band images as RGB and labeling the axes with proper geo-referenced extents.

from rasterio.plot import show
import rasterio
import geopandas as gpd

df = gpd.read_file(r"/home/bera/GIS/data/Indexrutor/index_10km.shp", 
                   bbox=(625263, 7323514, 677386, 7383838))
ax = df.boundary.plot(color="red", figsize=(15,15), linewidth=5)

demfile = r"/home/bera/GIS/data/nh_73_6.tif"
src = rasterio.open(demfile)
show(src, cmap="gray", ax=ax)

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.