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

Taras
  • 32,823
  • 4
  • 66
  • 137
BallpenMan
  • 1,217
  • 8
  • 19

1 Answers1

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.

Urban87
  • 616
  • 3
  • 11