USGS DEM Files: How to Load, Merge, and Crop with Python

QUICK SUCCESS DATA SCIENCE

A Digital Elevation Model (DEM) is a 3D digital representation of the earth’s surface. It records the height above sea level for various points, which may be measured using traditional surveying techniques, LIDAR, satellite imagery, aerial photography, or GPS measurements. DEMs are usually stored in raster format, where each pixel has an elevation value.

DEMs are essential tools in fields that use the physical landscape for planning and decision-making. Example uses include topographic mapping, hydrological modeling, urban planning, environmental studies, and creating realistic terrain models for virtual reality, gaming, and simulations.

The United States Geological Survey (USGS) provides free DEMs downloadable from the National Map Downloader. These are seamless raster files called DEM TIFFs or GeoTIFFs. While similar to standard photographic TIFFs, these files contain embedded Geospatial metadata.

In this Quick Success Data Science project, we’ll prepare a DEM TIFF for a selected geographical area. This won’t be as easy as downloading a file from a website. The existing USGS files will rarely coincide with the area of your study. And because they are memory intensive, you’ll want to trim them as much as possible. We’ll use the Rasterio library for this task and merge and crop the files to an area of interest (AOI) polygon.

---

Defining the Area of Interest

The DEM file we’ll build here will be used to study the watershed and stream networks of the Bayou Pierre, a river in southwest Mississippi. Thus, our AOI should encompass the river’s drainage basin, sometimes called the catchment area.

We’ll use the government map below to determine the size of this AOI. This will ensure our bounding box is big enough to encompass the river and the termination points of all its tributaries.

---

Downloading the DEM Data

USGS DEM TIFFs are large; the ones we’ll use here are around 60 MB each. We’ll use the National Map Downloader, previously mentioned, to download the files. Click here to start.

From the launch window select the 1 arc-second DEM under Elevation Products (3DEP), __ as shown in the following figure:

NOTE: This lower-resolution dataset will reduce memory requirements and speed processing.

Next, click the Enter Coords button and enter the following AOI coordinates:

Click the Add to Map button. You should see the following polygon appear:

This box should be large enough to encompass the river and its tributaries.

To fetch the DEM TIFFs that cover this box, click the blue Search Products button. You’ll see the results under the Products tab:

Four separate DEM TIFFs cover our AOI. To see the extent of each, hover your mouse over one of the items in the list. The bottom item yields the following map:

To create a DEM TIFF clipped to our AOI, we’ll need to download all four, merge them into one, and then crop the merged file to our bounding box.

To download the files, add them to the cart (don’t worry, they’re free) and then click the Cart tab. Click the four blue links under the Download column to complete the process.

Make a note of the download location, as you’ll need to add the path to the code.

---

The Code

With the DEM files downloaded, it’s time to merge and crop them with the Rasterio library. The following code, written in JupyterLab, is derived and modified from the Automating GIS Processes course (license info here).

The Rasterio Library

Rasterio is designed for reading, writing, and analyzing geospatial raster __ data. Raster data is a digital structure for capturing, storing, and representing spatial information. Values, representing data points on the ground, are stored in a matrix of cells or pixels. Example datasets are satellite images and DEMs.

To install Rasterio with pip use:

pip install rasterio

To install with conda use:

conda install -c conda-forge rasterio

Importing Libraries

You’ll also need to install Matplotlib (for plotting), GeoPandas (for working with shapefiles and polygons), and Fiona (for transforming and converting geospatial data between different formats). Installation instructions can be found in the previous hyperlinks. Fiona should install with GeoPandas.

Here are the imports:

import json
import matplotlib.pyplot as plt
import geopandas as gpd
from fiona.crs import from_epsg
from shapely.geometry import box

import rasterio
from rasterio.merge import merge
from rasterio.mask import mask
from rasterio.plot import show

Loading the DEM TIFFs and Checking the Metadata

The following code loads the USGS DEM TIFFs we downloaded previously and checks some of their stats. You’ll need to provide the proper path for your system.

# Load the DEM TIFFs and Check the stats:

# List of DEM TIFFs:
tifs = ['D:Pictures_on_Ddrainage DEMsUSGS_1_n32w091_20220601.tif', 
        'D:Pictures_on_Ddrainage DEMsUSGS_1_n32w092_20230609.tif',
        'D:Pictures_on_Ddrainage DEMsUSGS_1_n33w091_20221121.tif',
        'D:Pictures_on_Ddrainage DEMsUSGS_1_n33w092_20221121.tif']

# Loop through TIFFs and print stats:
for tif in tifs:
    with rasterio.open(tif) as src: 
        # Get the array shape and size:
        dem_array = src.read(1)  # Read the first band
        array_shape = dem_array.shape
        array_size = dem_array.size

        # Get the Coordinate Reference System (CRS):
        dem_crs = src.crs

    # Print the results
    print(f"DEM File: {tif}")
    print(f"Array Shape: {array_shape}")
    print(f"Array Size: {array_size}")
    print(f"CRS: {dem_crs}")
    print()

DEM File: D:Pictures_on_Ddrainage DEMsUSGS_1_n32w091_20220601.tif
Array Shape: (3612, 3612)
Array Size: 13046544
CRS: EPSG:4269

DEM File: D:Pictures_on_Ddrainage DEMsUSGS_1_n32w092_20230609.tif
Array Shape: (3612, 3612)
Array Size: 13046544
CRS: EPSG:4269

DEM File: D:Pictures_on_Ddrainage DEMsUSGS_1_n33w091_20221121.tif
Array Shape: (3612, 3612)
Array Size: 13046544
CRS: EPSG:4269

DEM File: D:Pictures_on_Ddrainage DEMsUSGS_1_n33w092_20221121.tif
Array Shape: (3612, 3612)
Array Size: 13046544
CRS: EPSG:4269

The files are the same size and shape and have the same coordinate reference system (CRS) for projecting onto a flat map.

Merging the DEM TIFFs

Now we use Rasterio to merge the four files into one. The code opens each TIFF, appends it to a list, and then passes the list to the Rasterio merge() method.

# Merge the individual DEM TIFFs into a single file:

# Assign list to hold opened tifs:
src_files_to_mosaic = []

# Create mosaic:
for tif in tifs:
    src = rasterio.open(tif)
    src_files_to_mosaic.append(src)

mosaic, out_trans = merge(src_files_to_mosaic)

The merge() method combines all the rasters, handling overlaps and stitching them together while ensuring they align correctly based on their geographic coordinates. It then returns two objects.

The mosaic variable represents the combined raster data. The out_trans variable is the transformation matrix for the mosaic, which includes geographical information like the origin, pixel size, and rotation.

This affine transformation matrix is crucial for maintaining the spatial integrity of the cropped DEM. It describes how to map pixel coordinates to geographic coordinates, adjusting the CRS to match the new spatial extent after cropping. This means that the geographic coordinates of the cropped DEM match the original DEM’s coordinate system.

Next, we copy the metadata from one of the TIFFs and embed it into a saved version of the combined TIFFs. We then show the merged image.

# Copy metadata from one of the DEM TIFFs
# (assumes all have the same metadata):

meta = src_files_to_mosaic[0].meta.copy()
meta.update({'driver': 'GTiff',
             'height': mosaic.shape[1],
             'width': mosaic.shape[2],
             'transform': out_trans})

# Save the mosaic as a new Geotiff file:
output_file = 'merged_mosaic.tif'
with Rasterio.open(output_file, "w", **meta) as dest:
    dest.write(mosaic)

# Show the mosaic
show(mosaic, cmap='terrain');

The first line copies the metadata from the first raster file in thesrc_files_to_mosaic list. Copying the metadata ensures that any changes made to meta do not affect the original metadata.

The next line updates the copied metadata with the file format (set to 'GTiff' (GeoTIFF)), then sets the number of rows (height) and columns (width) of the mosaic and updates the transformation matrix.

To save the data as a new file named _mergedmosaic.tif, use Rasterio to open the new file with the specified metadata (**meta). The 'w' indicates the file is opened in write mode.

We finish by showing the mosaic with the terrain colormap:

Cropping the Merged Files

We’ll use a rectangular polygon to crop the merged file to the extent of the Bayou Pierre watershed. We’ll load this box as a Geopandas’ GeoDataFrame and give it the same CRS as the DEM TIFFs. The following code will extract the box’s four coordinates in a Rasterio-compatible format.

# Prepare the bounding box for the AOI:

# Coordinates (lat-lon) for cropping:
minx, miny = -91.255, 31.6
maxx, maxy = -90.32, 32.25
bbox = box(minx, miny, maxx, maxy)

# Load bounding box into a GeoDataFrame:
geo = gpd.GeoDataFrame({'geometry': bbox}, 
                       index=[0], 
                       crs=from_epsg(4269))

def get_features(gdf):
    """Parse GDF features for use with Rasterio."""
    return [json.loads(gdf.to_json())['features'][0]['geometry']]

coords = get_features(geo)
print(coords)

Here’s the Rasterio-compatible output:

[{'type': 'Polygon', 'coordinates': [[[-90.32, 31.6], [-90.32, 32.25], [-91.255, 32.25], [-91.255, 31.6], [-90.32, 31.6]]]}]

Now we need to crop the merged file to this box. We accomplish this with the Rasterio mask() method, which takes the merged file and the box coordinates and returns an image along with the transform metadata. We then repeat the process of making a new output file and showing the results.

# Crop the merged DEM TIFF file:

# Path to file:
dem_file = 'merged_mosaic.tif'

# Open the DEM file:
with rasterio.open(dem_file) as src:
    # Clip the DEM file to the bounding box coords:
    out_img, out_transform = mask(dataset=src, 
                                  shapes=coords, 
                                  crop=True)

    # Update the metadata with the new transform and shape
    out_meta = src.meta.copy()
    out_meta.update({'driver': 'GTiff',
                     'height': out_img.shape[1],
                     'width': out_img.shape[2],
                     'transform': out_transform})

    # Save the cropped DEM:
    cropped_output = 'cropped_dem.tif'
    with rasterio.open(cropped_output, 'w', **out_meta) as dest:
        dest.write(out_img)

# Show the cropped DEM: 
with rasterio.open(cropped_output) as cropped_src: 
    fig, ax = plt.subplots(figsize=(10, 10)) 
    img = rasterio.plot.show(cropped_src, 
                             ax=ax, 
                             cmap='terrain') 
    # Add a color bar:
    cbar = plt.colorbar(img.get_images()[0], 
                        ax=ax, 
                        orientation='vertical', 
                        fraction=0.03, 
                        pad=0.04) 
    cbar.set_label('Elevation (meters)', 
                   rotation=270, 
                   labelpad=15)
    ax.set_title('DEM of Bayou Pierre Watershed')
    ax.set_xlabel('Longitude')
    ax.set_ylabel('Latitude')

That’s the main task completed. We now have a DEM TIFF cropped to our area of interest. From this we can generate derivative products, such as a flow accumulation diagram for the Bayou Pierre River system:

Next, we’ll embellish our DEM TIFF with some reference lines and markers.

Drawing on the DEM TIFF

As a bonus project, let’s add the boundary for Claiborne County and mark the highest elevation within the county.

For the boundary polygon, we’ll use a shapefile. This is a widely used format for storing vector data in geographic information systems. You can read more about it here:

Shape Up Your Maps with Shapefiles

Below is the complete annotated code for this operation. You can download the county shapefile (_cb_2022_us_county500k.zip) from the US Census Bureau.

# Plot a DEM with a county boundary.
# Mark the highest point in the county.

import rasterio
import rasterio.mask
import geopandas as gpd
import matplotlib.pyplot as plt
from rasterio.plot import show
from shapely.geometry import mapping, box
import numpy as np
from matplotlib.lines import Line2D

# Load the county shapefile:
COUNTY_NAME = 'Claiborne'
SHAPEFILE_NAME = 'cb_2022_us_county_500k.zip'

county_gdf = gpd.read_file(SHAPEFILE_NAME)
county_gdf.head()

# Examine the GDF and use the proper column name for county name:
# county_gdf.head()
county_col = 'NAME'  # This name is unique to shapefile!

county = county_gdf[(county_gdf[county_col] == COUNTY_NAME)]

# Read the DEM TIFF file:
with rasterio.open('cropped_dem.tif') as cropped_src:

    # Ensure the shapefile is in the same CRS as the DEM:
    if county.crs != cropped_src.crs:
        county = county.to_crs(cropped_src.crs)

    # Get the bounding box of the DEM:
    dem_bounds = cropped_src.bounds
    dem_bbox = box(dem_bounds.left, 
                   dem_bounds.bottom, 
                   dem_bounds.right, 
                   dem_bounds.top)

    # Clip the shapefile to the bounding box of the DEM.
    # In case the county is larger than the DEM coverage.
    county_clipped = county.clip(dem_bbox)

    # Mask the DEM using the polygon:
    geoms = [mapping(geom) for geom in county_clipped.geometry]
    out_image, out_transform = rasterio.mask.mask(cropped_src, 
                                                  geoms, 
                                                  crop=True)

    # Extract the max elevation value and its location:
    max_elev = np.max(out_image)
    max_index = np.argmax(out_image)
    max_coords = np.unravel_index(max_index, out_image.shape[1:])

    # Print the max elevation for reference:
    print(f"Maximum elevation in county is: {max_elev:.2f} m")

    # Convert array coordinates to geospatial coordinates
    max_lon, max_lat = rasterio.transform.xy(out_transform, 
                                             max_coords[0], 
                                             max_coords[1])

    # Plot the DEM and the county polygon
    fig, ax = plt.subplots(figsize=(10, 10))
    img = show(cropped_src, ax=ax, cmap='terrain')

    # Create a colorbar with a legend
    cbar = plt.colorbar(img.get_images()[0], 
                        ax=ax, 
                        orientation='vertical', 
                        fraction=0.03, 
                        pad=0.04)
    cbar.set_label('Elevation (meters)', rotation=270, labelpad=15)

    # Plot the county polygon
    county_clipped.plot(ax=ax, 
                        facecolor='none', 
                        edgecolor='red', 
                        linewidth=1)

    # Plot the maximum elevation point
    ax.plot(max_lon, max_lat, 
            'bo', 
            markersize=10, 
            label='Max Elevation')

    # Create custom legend handles
    custom_lines = [Line2D([0], [0], 
                           color='red', 
                           lw=2, 
                           label=f'{COUNTY_NAME} County Boundary'),
                    Line2D([0], [0], 
                           marker='o', 
                           color='w', 
                           markerfacecolor='blue', 
                           markersize=10, 
                           label='Max Elevation')]

    # Set title and labels
    plt.title(f'USGS DEM with {COUNTY_NAME} County Boundary')
    plt.xlabel('Longitude')
    plt.ylabel('Latitude')

    # Show the plot with the legend
    plt.legend(handles=custom_lines)
    plt.show()

The highest point in Claiborne County is 135 m (443 ft). Growing up, I had always heard that the highest point was in the Sunset Hill area, circled in blue in the figure below. Sunset Hill, however, is only 119 m (390′). This urban legend probably arose because Sunset Hill is close to the county seat, meaning more people knew about it. The true highest point is way off in the boondocks!

Besides being fun and useful, DEMs also make for some great art. I once saw one etched in a block of thick glass, which was quite fetching.

---

Summary

Digital Elevation Models (DEMs) are digital files that capture the height above sea level for a part of the earth’s surface. The US Geological Survey (USGS) provides these files for free for areas within the country. They come as DEM TIFFs (GeoTIFFS) in raster format.

DEMs are useful for urban planning, floodplain mapping, environmental studies, and more. They can also be used to model terrains in virtual reality applications.

In this project, we learned how to download multiple DEM files, merge them into a single mosaic, and crop the mosaic with a boundary polygon. We also learned how to query the file for elevation data and annotate it with a political boundary.

---

Data Permissions

[1] Bayou Pierre Sub-basin map: Permission to use granted by: Stephen Champlin, RPG, Geospatial Resources Division/Flood Mapping Director, Office of Geology, Mississippi Department of Environmental Quality (MDEQ), http://geology.deq.ms.gov/floodmaps.

---

Thanks!

Thanks for reading and please follow me for more Quick Success Data Science projects in the future. If you found this article useful, please clap (up to 50 claps are allowed), highlight text, or leave a comment. Your engagement helps authors earn more, and we greatly appreciate it.