How to Get Surface Elevation for a Specific Latitude and Longitude With Python

Overview

A simple way to get surface elevation for a latitude/longitude point is to sample a Digital Elevation Model. Here we use the Copernicus DEM GLO-30, a global 30 m Digital Surface Model available as Cloud Optimized GeoTIFFs on AWS. It represents the Earth’s surface, including terrain, buildings, and vegetation; elevation values are in meters. (Open Data on AWS)

Install packages

1
conda install -c conda-forge rasterio numpy

Example point

In this example, we retrieve the elevation around Mont Blanc, located in the French Alps near the border between France and Italy.

1
2
3
4
5
6
7
8
import os
import rasterio
import numpy as np

os.environ["AWS_NO_SIGN_REQUEST"] = "YES"

lon = 6.8652
lat = 45.8326

Open the correct Copernicus DEM tile

Copernicus DEM tiles are organized by 1° × 1° folders on AWS.

For latitude 45.8326 and longitude 6.8652, the tile is:

1
2
3
4
5
6
7
tile = "Copernicus_DSM_COG_10_N45_00_E006_00_DEM"

url = (
    f"/vsis3/copernicus-dem-30m/"
    f"{tile}/"
    f"{tile}.tif"
)

Open the tile:

1
2
with rasterio.open(url) as src:
    print(src.profile)

Sample elevation at one point

1
2
3
4
with rasterio.open(url) as src:
    elevation_m = list(src.sample([(lon, lat)]))[0][0]

print(elevation_m)

Example output:

1
4797.0034

This means:

1
2
4797.0034 meters above sea level
 4.78 km

Rasterio can open local files, URLs, or GDAL-supported virtual paths such as /vsis3/, and sample() extracts raster values at coordinate locations. (rasterio.readthedocs.io)

Visualize point on a satellite map with Folium

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
import folium

# Mont Blanc coordinates
lon = 6.8652
lat = 45.8326

# Create map
m = folium.Map(
    location=[lat, lon],
    zoom_start=12,
    tiles=None
)

# Add Esri satellite imagery
folium.TileLayer(
    tiles=(
        "https://server.arcgisonline.com/ArcGIS/rest/services/"
        "World_Imagery/MapServer/tile/{z}/{y}/{x}"
    ),
    attr="Esri World Imagery",
    name="Esri Satellite",
    overlay=False,
    control=True,
).add_to(m)

# Add marker
folium.CircleMarker(
    location=[lat, lon],
    radius=8,
    color="red",
    fill=True,
    fill_color="red",
    fill_opacity=0.9,
    popup=(
        f"<b>Mont Blanc</b><br>"
        f"Latitude: {lat}<br>"
        f"Longitude: {lon}"
    ),
).add_to(m)

# Optional layer control
folium.LayerControl().add_to(m)

# Display map
m

How to Get Surface Elevation for a Specific Latitude and Longitude With Python
How to Get Surface Elevation for a Specific Latitude and Longitude With Python

Sample multiple points

1
2
3
4
5
6
7
8
9
lons = np.array([-104.98, -104.50])
lats = np.array([40.20, 40.80])

coords = list(zip(lons, lats))

with rasterio.open(url) as src:
    elevations_m = np.array([v[0] for v in src.sample(coords)])

print(elevations_m)

Example output:

1
[1515.6829 1583.732 ]

Add elevation to a DataFrame

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
import pandas as pd

df = pd.DataFrame({
    "longitude": [-104.98, -104.50],
    "latitude": [40.20, 40.80],
})

coords = list(zip(df["longitude"], df["latitude"]))

with rasterio.open(url) as src:
    df["elevation_m"] = [v[0] for v in src.sample(coords)]

df["elevation_km"] = df["elevation_m"] / 1000

print(df)

Function for one tile

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
def sample_elevation_from_tile(lons, lats, tile):
    url = (
        f"/vsis3/copernicus-dem-30m/"
        f"{tile}/"
        f"{tile}.tif"
    )

    coords = list(zip(lons, lats))

    with rasterio.open(url) as src:
        elevations = np.array([v[0] for v in src.sample(coords)])

    return elevations

Usage:

1
2
3
4
5
6
7
tile = "Copernicus_DSM_COG_10_N40_00_W105_00_DEM"

df["elevation_m"] = sample_elevation_from_tile(
    df["longitude"].values,
    df["latitude"].values,
    tile
)

Important note

This example works only when all points fall inside the same DEM tile. For larger datasets, you need to group points by tile name first, then sample each tile separately.

Visualize points on a satellite map with Folium

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
import folium

# Center map
center_lat = lats.mean()
center_lon = lons.mean()

m = folium.Map(
    location=[center_lat, center_lon],
    zoom_start=12,
    tiles=None
)

# Satellite imagery
folium.TileLayer(
    tiles=(
        "https://server.arcgisonline.com/ArcGIS/rest/services/"
        "World_Imagery/MapServer/tile/{z}/{y}/{x}"
    ),
    attr="Esri World Imagery",
    name="Esri Satellite",
    overlay=False,
    control=True,
).add_to(m)

# Add points
for i, (lat, lon, elev) in enumerate(
    zip(lats, lons, elevations_m)
):

    popup_text = (
        f"<b>Point {i}</b><br>"
        f"Latitude: {lat:.4f}<br>"
        f"Longitude: {lon:.4f}<br>"
        f"Elevation: {elev:.1f} m"
    )

    folium.CircleMarker(
        location=[lat, lon],
        radius=7,
        color="red",
        fill=True,
        fill_color="red",
        fill_opacity=0.9,
        popup=popup_text,
    ).add_to(m)

folium.LayerControl().add_to(m)

m

Save map as HTML

1
m.save("mont_blanc_elevation_map.html")

Find the Correct Copernicus DEM Tile for a Latitude and Longitude

The Copernicus DEM is organized into 1° × 1° geographic tiles.

Each tile name follows this structure:

1
Copernicus_DSM_COG_10_<LAT>_00_<LON>_00_DEM

Examples:

1
2
Copernicus_DSM_COG_10_N45_00_E006_00_DEM
Copernicus_DSM_COG_10_N40_00_W105_00_DEM

Step 1 — Determine tile indices

For a coordinate:

1
2
Latitude  = 45.8326
Longitude = 6.8652

take the integer floor:

1
2
3
4
5
6
7
8
9
import numpy as np

lat = 45.8326
lon = 6.8652

lat_floor = int(np.floor(lat))
lon_floor = int(np.floor(lon))

print(lat_floor, lon_floor)

Output:

1
45 6

Step 2 — Convert to Copernicus naming format

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
if lat_floor >= 0:
    lat_str = f"N{lat_floor:02d}"
else:
    lat_str = f"S{abs(lat_floor):02d}"

if lon_floor >= 0:
    lon_str = f"E{lon_floor:03d}"
else:
    lon_str = f"W{abs(lon_floor):03d}"

print(lat_str, lon_str)

Output:

1
N45 E006

Step 3 — Build the tile name

1
2
3
4
5
6
7
tile = (
    f"Copernicus_DSM_COG_10_"
    f"{lat_str}_00_"
    f"{lon_str}_00_DEM"
)

print(tile)

Output:

1
Copernicus_DSM_COG_10_N45_00_E006_00_DEM

Complete function

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import numpy as np

def get_copernicus_tile(lat, lon):

    lat_floor = int(np.floor(lat))
    lon_floor = int(np.floor(lon))

    # Latitude string
    if lat_floor >= 0:
        lat_str = f"N{lat_floor:02d}"
    else:
        lat_str = f"S{abs(lat_floor):02d}"

    # Longitude string
    if lon_floor >= 0:
        lon_str = f"E{lon_floor:03d}"
    else:
        lon_str = f"W{abs(lon_floor):03d}"

    tile = (
        f"Copernicus_DSM_COG_10_"
        f"{lat_str}_00_"
        f"{lon_str}_00_DEM"
    )

    return tile

Example: Mont Blanc

1
2
3
4
5
6
lat = 45.8326
lon = 6.8652

tile = get_copernicus_tile(lat, lon)

print(tile)

Output:

1
Copernicus_DSM_COG_10_N45_00_E006_00_DEM

Example: Colorado

1
2
3
4
5
6
lat = 40.20
lon = -104.98

tile = get_copernicus_tile(lat, lon)

print(tile)

Output:

1
Copernicus_DSM_COG_10_N40_00_W105_00_DEM

Build the AWS URL automatically

1
2
3
4
5
6
7
8
9
tile = get_copernicus_tile(lat, lon)

url = (
    f"/vsis3/copernicus-dem-30m/"
    f"{tile}/"
    f"{tile}.tif"
)

print(url)

Output:

1
2
3
/vsis3/copernicus-dem-30m/
Copernicus_DSM_COG_10_N45_00_E006_00_DEM/
Copernicus_DSM_COG_10_N45_00_E006_00_DEM.tif

Each DEM tile:

  • spans 1 degree latitude
  • spans 1 degree longitude
  • contains a 30 m resolution elevation grid

Important note about negative coordinates

Coordinate Prefix
Positive latitude N
Negative latitude S
Positive longitude E
Negative longitude W

Examples:

Location Tile
France N45 E006
Colorado N40 W105
Chile S33 W071

References

Resource Why useful
Copernicus DEM on AWS Dataset access and S3 structure
Copernicus DEM Product Handbook Vertical reference and product details
Rasterio documentation Reading GeoTIFFs and sampling raster values
Image

of