This guide has 2 routes. Pick one, or read both. Each route is complete on its own.
A digital elevation model, or DEM, is a raster in which each cell holds a height. It is the starting point for slope, aspect, hillshade, watersheds, viewsheds and flood work. Free global DEMs exist at 30 m, and the best known are SRTM, ALOS World 3D and the Copernicus DEM. This guide uses the Copernicus DEM GLO-30, which covers the world at 30 m and sits on open servers, so you can fetch a tile with a link and no account.
Know what you are downloading
DEM is the general word for any raster of heights. The two common kinds differ in what the height is measured on.
A DSM, a digital surface model, records the top of whatever the satellite sees first: the roof of a house, the canopy of a forest, the ground in a bare field. A DTM, a digital terrain model, records the bare ground, with buildings and trees removed. The Copernicus DEM is a DSM. In a forest it shows the canopy, and in a city it shows the roofs. The study note DEM, DSM and DTM explains the difference and why it matters for flood and slope work.
The Copernicus DEM comes in two resolutions. GLO-30 has a cell of 1 arc second, about 30 m. GLO-90 has 3 arc seconds, about 90 m. Use GLO-30 unless the area is huge. Each tile covers 1 degree by 1 degree of latitude and longitude. A small number of GLO-30 tiles, for a few countries, are not released to the public, and the ocean has no tiles at all. Where tiles are missing, GLO-90 may still exist.
Read a tile name
A tile name looks like Copernicus_DSM_COG_10_N30_00_E075_00_DEM. It reads:
Copernicus_DSM_COGis the product, a surface model stored as a Cloud Optimized GeoTIFF.10is the cell size in arc seconds times ten, so10is GLO-30. GLO-90 tiles use30.N30_00is the latitude of the south-west corner of the tile, 30 degrees north.E075_00is the longitude of the south-west corner, 75 degrees east.
The tile therefore runs from 30 N to 31 N and 75 E to 76 E. Jagraon in Punjab, at about 30.79 N and 75.47 E, is inside it. To find the tile for your place, round your latitude and longitude down to whole degrees. Use N for north and S for south, E for east and W for west.
Route A: download a tile with a link
The tiles sit in a public bucket on Amazon's servers. A link to the tile file downloads it, and no account is needed. Build the link from the tile name:
https://copernicus-dem-30m.s3.amazonaws.com/<tile name>/<tile name>.tif
For the tile above, paste this into your browser:
https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N30_00_E075_00_DEM/Copernicus_DSM_COG_10_N30_00_E075_00_DEM.tif
The file is about 40 MB. For GLO-90, replace copernicus-dem-30m with copernicus-dem-90m and 10 in the name with 30. Open the downloaded file in QGIS. It is a GeoTIFF in latitude and longitude (EPSG:4326), and the heights are in metres.
If you prefer a portal, two serve global DEMs for an area that you choose.
OpenTopography has a page for the Copernicus GLO-30 Digital Elevation Model. Its overview says in so many words that the product is a Digital Surface Model. The page has a map headed Select area of data to process, and four numbered sections. 1. Coordinates shows the coordinate system, WGS84 (EPSG:4326) horizontally and WGS84 Ellipsoid with the EGM2008 geoid vertically. It has boxes for the minimum and maximum X and Y, which you can type if you tick Manually enter selection coordinates, and a link to Validate coordinates and estimate area. 2. Data Output Formats offers Arc ASCII Grid, GeoTiff and IMG. 3. Layer Types has a checkbox for the Digital Surface Model (DSM). 4. Raster Visualization offers a hillshade and a colour map. A Submit button sends the request. The page greets you as a guest, with a Sign In link, and the download step may ask you to sign in. That step was not tested.
USGS EarthExplorer has a Digital Elevation category, which holds SRTM data sets such as Shuttle Radar Topography Mission 1 Arc-Second Digital Terrain Elevation Data - Global and SRTM Void Filled. Searching needs no account, and downloading needs a free one. Steps 1 to 3 of the route in the Landsat guide apply, with this category in place of Landsat. NASA Earthdata also serves SRTM and needs an account.
The link above is faster when you know the tile.
Route B: Python
The script below finds the tile for a box, reads only your area from the online file, checks the heights, and saves the result. You need Python with requests, numpy and rasterio.
import requests
# Step 1: describe your area. The box is west, south, east, north in longitude and latitude.
study_area = [75.40, 30.74, 75.55, 30.84] # a box of about 14 km by 11 km near Jagraon, Punjab
# Step 2: ask the Earth Search catalogue which Copernicus DEM tiles touch the box
search = {"collections": ["cop-dem-glo-30"], "bbox": study_area, "limit": 10}
url = "https://earth-search.aws.element84.com/v1/search"
tiles = requests.post(url, json=search, timeout=60).json()["features"]
# Step 3: print each tile name and the address of its file
for tile in tiles:
print(tile["id"])
print(" ", tile["assets"]["data"]["href"])
Copernicus_DSM_COG_10_N30_00_E075_00_DEM
s3://copernicus-dem-30m/Copernicus_DSM_COG_10_N30_00_E075_00_DEM/Copernicus_DSM_COG_10_N30_00_E075_00_DEM.tif
The catalogue gives the address in an Amazon-specific form that begins with s3://. A normal web address for the same file replaces the start with https://<bucket>.s3.amazonaws.com/. The next step does that and reads your area.
import numpy as np
import rasterio
from rasterio.windows import from_bounds
# Step 1: turn the s3:// address into a web address that rasterio can open
def web_address(s3_address):
bucket, key = s3_address.replace("s3://", "").split("/", 1)
return f"https://{bucket}.s3.amazonaws.com/{key}"
# Step 2: open the tile online and read only the window around your box. The file is in degrees, so the box needs no conversion.
address = web_address(tiles[0]["assets"]["data"]["href"])
with rasterio.open(address) as src:
window = from_bounds(*study_area, src.transform).round_offsets().round_lengths()
heights = src.read(1, window=window).astype("float32")
transform = src.window_transform(window)
crs = src.crs
print("tile CRS: EPSG", crs.to_epsg(), "| cell size:", round(src.res[0], 6), "degrees | nodata:", src.nodata)
# Step 3: the facts about what you read
print("window:", heights.shape[1], "x", heights.shape[0], "cells")
print("height: lowest", round(float(heights.min()), 1), "m | highest", round(float(heights.max()), 1), "m | mean", round(float(heights.mean()), 1), "m")
tile CRS: EPSG 4326 | cell size: 0.000278 degrees | nodata: None
window: 540 x 360 cells
height: lowest 227.5 m | highest 253.2 m | mean 233.6 m
The cell is about 0.000278 degrees, which is one arc second. In degrees the cell is a constant, but on the ground it is about 30 m from north to south and a little less from east to west, because meridians come together towards the poles. This matters for the next step, and the "Common mistakes" section returns to it.
Check the heights
A quick check catches most errors. Jagraon is a flat town on the Punjab plain at roughly 230 to 250 m above sea level. The window should show something like that, with a spread of a few tens of metres at most.
# Step 1: the spread of heights in the window. A flat plain has a small spread.
spread = float(heights.max() - heights.min())
print("spread:", round(spread, 1), "m")
# Step 2: a rough slope check. Height change per cell in metres, divided by the size of a cell on the ground.
cell_north_south = 30.9 # one arc second of latitude is about 30.9 m
cell_east_west = 30.9 * np.cos(np.radians(30.79)) # east-west cells are narrower, by the cosine of the latitude
dz_dy, dz_dx = np.gradient(heights, cell_north_south, cell_east_west)
slope = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))
# Step 3: the typical and the steepest slope
print("median slope:", round(float(np.median(slope)), 2), "degrees | steepest:", round(float(slope.max()), 1), "degrees")
spread: 25.7 m
median slope: 0.37 degrees | steepest: 20.5 degrees
On a plain the median slope is well under one degree. The steepest cells, here over 20 degrees, sit at the edges of buildings, canal banks and embankments. That is the surface model at work: it records the top of a house and the ground beside it in neighbouring cells. If the median slope is tens of degrees, or the steepest is near 90, the height and the cell size are in different units, which is the most common DEM mistake.
Join two tiles
If your box crosses a tile edge, the catalogue returns two tiles, and one read gives only half your area. rasterio.merge reads and joins tiles for exactly the box you give it. This box straddles the 75 E line between two tiles.
from rasterio.merge import merge
# Step 1: a box that crosses a tile edge at 75 degrees east
edge_box = [74.90, 30.85, 75.10, 30.95]
search = {"collections": ["cop-dem-glo-30"], "bbox": edge_box, "limit": 10}
edge_tiles = requests.post(url, json=search, timeout=60).json()["features"]
print(len(edge_tiles), "tiles:", [t["id"].replace("Copernicus_DSM_COG_10_", "").replace("_DEM", "") for t in edge_tiles])
# Step 2: open them all online and let rasterio merge them into the box
files = [rasterio.open(web_address(t["assets"]["data"]["href"])) for t in edge_tiles]
mosaic, mosaic_transform = merge(files, bounds=tuple(edge_box))
for f in files:
f.close()
# Step 3: one array covering the whole box
print("joined window:", mosaic.shape[2], "x", mosaic.shape[1], "cells | mean height:", round(float(mosaic.mean()), 1), "m")
2 tiles: ['N30_00_E075_00', 'N30_00_E074_00']
joined window: 720 x 360 cells | mean height: 212.8 m
Save the result
# Step 1: write the heights to a GeoTIFF that keeps the position and CRS of the window
with rasterio.open(
"dem_cop30.tif", "w", driver="GTiff",
height=heights.shape[0], width=heights.shape[1], count=1,
dtype="float32", crs=crs, transform=transform, nodata=-9999, compress="deflate",
) as dst:
dst.write(heights, 1)
# Step 2: read it back to check
with rasterio.open("dem_cop30.tif") as check:
print("saved:", check.width, "x", check.height, "cells | CRS EPSG", check.crs.to_epsg(), "| data type", check.dtypes[0])
saved: 540 x 360 cells | CRS EPSG 4326 | data type float32
Vertical datum and CRS
Two details decide whether the DEM fits with your other data.
The horizontal CRS of the tile is geographic, EPSG:4326, so the cell size is in degrees. For slope, aspect, hillshade and area, reproject the DEM to a projected CRS in metres, such as the UTM zone of your place, first. Reproject before you calculate and not after, otherwise the slope has the problem the check above warned about.
The heights are above a geoid, a model of mean sea level. The Copernicus DEM uses the EGM2008 geoid, so heights are metres above mean sea level and not above the WGS 84 ellipsoid. GPS devices often report ellipsoid heights, which differ from geoid heights by tens of metres. Do not subtract one from the other without converting.
Common mistakes
- Treating a DSM as bare ground. In a forest or a city the heights include the canopy and the roofs.
- Calculating slope on a DEM in degrees with metres as the height unit. The slope comes out near 90 degrees.
- Reading one tile for a box that crosses a tile edge, and getting half the area.
- Taking a tile for the wrong corner of the world, for example
E075when the place is at 75 W. - Mixing GPS ellipsoid heights with DEM geoid heights.
- Downloading a country's worth of tiles when you need a district. Read a window instead.