17. Working with Raster Data Using Rasterio#
17.1. Introduction#
17.1.1. What is Raster Data?#
17.1.2. Why Use Rasterio?#
17.2. Learning Objectives#
17.3. Installing Rasterio#
# %pip install rasterio pygis
import rasterio
import rasterio.plot
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
17.4. Reading Raster Data#
17.4.1. Opening Raster Files#
raster_path = (
"https://github.com/opengeos/datasets/releases/download/raster/dem_90m.tif"
)
src = rasterio.open(raster_path)
print(src)
17.4.2. Understanding Raster Metadata#
17.4.2.1. Basic File Information#
print(f"File name: {src.name}")
print(f"File mode: {src.mode}")
print("Raster metadata:")
for key, value in src.meta.items():
print(f"{key}: {value}")
17.4.2.2. Spatial Properties#
print(f"Coordinate Reference System: {src.crs}")
print(f"Pixel size (x, y): {src.res}")
print(f"Raster dimensions: {src.width} x {src.height} pixels")
print(f"Geographic bounds: {src.bounds}")
17.4.2.3. Data Properties#
print(f"Data types: {src.dtypes}")
print(f"Number of bands: {src.count}")
17.4.3. The Affine Transform#
print("Affine transform:")
print(src.transform)
17.5. Visualizing Raster Data#
17.5.1. Basic Raster Visualization#
rasterio.plot.show(src)
17.5.2. Understanding Color Maps#
fig, ax = plt.subplots(figsize=(8, 8))
rasterio.plot.show(src, cmap="terrain", ax=ax, title="Digital Elevation Model (DEM)")
plt.show()
17.5.3. Adding Colorbars#
elev_band = src.read(1)
plt.figure(figsize=(8, 8))
plt.imshow(elev_band, cmap="terrain")
plt.colorbar(label="Elevation (meters)", shrink=0.5)
plt.title("DEM with Terrain Colormap")
plt.show()
17.5.4. Visualizing Multiple Bands#
raster_path = "https://github.com/opengeos/datasets/releases/download/raster/LC09_039035_20240708_90m.tif"
src = rasterio.open(raster_path)
17.5.4.1. Single Band Visualization#
rasterio.plot.show((src, 1), cmap="gray", title="Band 1")
17.5.4.2. RGB Composite#
nir_band = src.read(5)
red_band = src.read(4)
green_band = src.read(3)
# Stack the bands into a single array
rgb = np.dstack((nir_band, red_band, green_band)).clip(0, 1)
# Plot the stacked array
plt.figure(figsize=(8, 8))
plt.imshow(rgb)
plt.title("Bands NIR, Red, and Green combined")
plt.show()
17.5.4.3. Creating a Multi-Panel Plot#
band_names = ["Coastal Aerosol", "Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"]
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(8, 10))
axes = axes.flatten() # Flatten the 2D array of axes to 1D for easy iteration
for band in range(1, src.count):
data = src.read(band)
ax = axes[band - 1]
im = ax.imshow(data, cmap="gray", vmin=0, vmax=0.5)
ax.set_title(f"Band {band_names[band - 1]}")
fig.colorbar(im, ax=ax, label="Reflectance", shrink=0.5)
plt.tight_layout()
plt.show()
17.5.5. Overlaying Vector Data#
# Load raster data
raster_path = (
"https://github.com/opengeos/datasets/releases/download/raster/dem_90m.tif"
)
src = rasterio.open(raster_path)
# Load vector data
vector_path = (
"https://github.com/opengeos/datasets/releases/download/places/dem_bounds.geojson"
)
gdf = gpd.read_file(vector_path)
gdf = gdf.to_crs(src.crs) # Ensure same CRS as raster
# Create the plot
fig, ax = plt.subplots(figsize=(8, 8))
rasterio.plot.show(src, cmap="terrain", ax=ax, title="DEM with Vector Overlay")
gdf.plot(ax=ax, edgecolor="red", facecolor="none", linewidth=2)
plt.show()
17.6. Accessing and Manipulating Raster Bands#
17.6.1. Stacking Multiple Bands#
raster_path = "https://github.com/opengeos/datasets/releases/download/raster/LC09_039035_20240708_90m.tif"
src = rasterio.open(raster_path)
nir_band = src.read(5)
red_band = src.read(4)
green_band = src.read(3)
# Stack the bands into a single array
rgb = np.dstack((nir_band, red_band, green_band)).clip(0, 1)
print(rgb.shape)
17.6.2. Basic Band Math (NDVI Calculation)#
# NDVI Calculation: NDVI = (NIR - Red) / (NIR + Red)
ndvi = (nir_band - red_band) / (nir_band + red_band)
ndvi = ndvi.clip(-1, 1)
plt.figure(figsize=(8, 8))
plt.imshow(ndvi, cmap="RdYlGn", vmin=-1, vmax=1)
plt.colorbar(label="NDVI", shrink=0.75)
plt.title("NDVI")
plt.xlabel("Column #")
plt.ylabel("Row #")
plt.show()
17.7. Writing Raster Data#
with rasterio.open(raster_path) as src:
profile = src.profile
print(profile)
profile.update(dtype=rasterio.float32, count=1, compress="lzw")
print(profile)
output_raster_path = "ndvi.tif"
with rasterio.open(output_raster_path, "w", **profile) as dst:
dst.write(ndvi, 1)
print(f"Raster data has been written to {output_raster_path}")
17.8. Clipping Raster Data#
17.8.1. Clipping with a Bounding Box#
raster_path = "https://github.com/opengeos/datasets/releases/download/raster/LC09_039035_20240708_90m.tif"
src = rasterio.open(raster_path)
data = src.read()
data.shape
subset = data[:, 900:1400, 700:1200].clip(0, 1)
rgb_subset = np.dstack((subset[4], subset[3], subset[2]))
rgb_subset.shape
# Plot the stacked array
plt.figure(figsize=(8, 8))
plt.imshow(rgb_subset)
plt.title("Las Vegas, NV")
plt.show()
from rasterio.windows import Window
from rasterio.transform import from_bounds
# Assuming subset and src are already defined
# Define the window of the subset (replace with actual window coordinates)
window = Window(col_off=700, row_off=900, width=500, height=500)
# Calculate the bounds of the window
window_bounds = rasterio.windows.bounds(window, src.transform)
# Calculate the new transform based on the window bounds
new_transform = from_bounds(*window_bounds, window.width, window.height)
with rasterio.open(
"las_vegas.tif",
"w",
driver="GTiff",
height=subset.shape[1],
width=subset.shape[2],
count=subset.shape[0],
dtype=subset.dtype,
crs=src.crs,
transform=new_transform,
compress="lzw",
) as dst:
dst.write(subset)
17.8.2. Clipping with a Vector Dataset#
import fiona
import rasterio.mask
geojson_path = "https://github.com/opengeos/datasets/releases/download/places/las_vegas_bounds_utm.geojson"
bounds = gpd.read_file(geojson_path)
fig, ax = plt.subplots()
rasterio.plot.show(src, ax=ax)
bounds.plot(ax=ax, edgecolor="red", facecolor="none")
with fiona.open(geojson_path, "r") as f:
shapes = [feature["geometry"] for feature in f]
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
out_meta = src.meta
out_meta.update(
{
"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform,
}
)
with rasterio.open("las_vegas_clip.tif", "w", **out_meta) as dst:
dst.write(out_image)