Cropping and Plotting Raster data in python using rasterio¶
This code uses the libraries shapely, rasterio, and matplotlib to extract a section of data from a raster image file (ETopo.tif), and plot it as a map of terrain elevations. You can find complete details on how to open and plot raster here.
Shapely is a library for manipulation and analysis of planar geometric objects, Rasterio is a library for reading and writing geospatial raster data, and Matplotlib is a plotting library for creating static, animated, and interactive visualizations in Python.
from shapely.geometry import MultiPolygon, Polygon
import rasterio.mask
import rasterio
import matplotlib.pyplot as plt
MultiPolygon and Polygon objects from shapely.geometry are used to define a bounding box (bbox) with a set of minimum and maximum values for the x and y coordinates of the area of interest.
xmin,xmax,ymin,ymax=60,100,0,40 # set the region
bbox=MultiPolygon([Polygon([[xmin, ymin], [xmin, ymax], [xmax, ymax], [xmax, ymin]])])
The rasterio.mask function is used to extract the data from the raster file (ETopo.tif) within the bounding box, using the src object created with rasterio.open.
# Open the raster data using rasterio
with rasterio.open("ETopo.tif") as src:
data,_=rasterio.mask.mask(src,shapes=bbox,crop=True)
matplotlib.pyplot is used to create a figure and axis objects, plot the extracted data as an image using the imshow method, and add a color bar to the plot with labels and extensions.
Finally, the show method is called to display the plot.
fig = plt.figure(figsize=[12,8])
# Plot the raster data using matplotlib
ax = fig.add_axes([0, 0, 1, 1])
raster_image=ax.imshow(data[0,:,:],cmap="terrain",vmax=5000,vmin=-4000)
fig.colorbar(raster_image, ax=ax,label="Elevation (in m) ",orientation='vertical',extend='both',shrink=0.5)
plt.show()
xmin,xmax,ymin,ymax=-120,-60,10,50
bbox=MultiPolygon([Polygon([[xmin, ymin], [xmin, ymax], [xmax, ymax], [xmax, ymin]])])
# Open the raster data using rasterio
with rasterio.open("ETopo.tif") as src:
data,_=rasterio.mask.mask(src,shapes=bbox,crop=True)
fig = plt.figure(figsize=[12,8])
# Plot the raster data using matplotlib
ax = fig.add_axes([0, 0, 1, 1])
raster_image=ax.imshow(data[0,:,:],cmap="terrain",vmax=5000,vmin=-3000)
fig.colorbar(raster_image, ax=ax,label="Elevation (in m) ",orientation='vertical',extend='both',shrink=0.5)
plt.show()
