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
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.
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:
tile = "Copernicus_DSM_COG_10_N45_00_E006_00_DEM"
url = (
f "/vsis3/copernicus-dem-30m/"
f " { tile } /"
f " { tile } .tif"
)
Open the tile:
with rasterio . open ( url ) as src :
print ( src . profile )
Sample elevation at one point
with rasterio . open ( url ) as src :
elevation_m = list ( src . sample ([( lon , lat )]))[ 0 ][ 0 ]
print ( elevation_m )
Example output:
This means:
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
Sample multiple points
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:
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:
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
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:
Copernicus_DSM_COG_10_ < LAT > _00_ < LON > _00_DEM
Examples:
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:
Latitude = 45.8326
Longitude = 6.8652
take the integer floor:
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:
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:
Step 3 — Build the tile name
tile = (
f "Copernicus_DSM_COG_10_"
f " { lat_str } _00_"
f " { lon_str } _00_DEM"
)
print ( tile )
Output:
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
lat = 45.8326
lon = 6.8652
tile = get_copernicus_tile ( lat , lon )
print ( tile )
Output:
Copernicus_DSM_COG_10_N45_00_E006_00_DEM
Example: Colorado
lat = 40.20
lon = - 104.98
tile = get_copernicus_tile ( lat , lon )
print ( tile )
Output:
Copernicus_DSM_COG_10_N40_00_W105_00_DEM
Build the AWS URL automatically
tile = get_copernicus_tile ( lat , lon )
url = (
f "/vsis3/copernicus-dem-30m/"
f " { tile } /"
f " { tile } .tif"
)
print ( url )
Output:
/ vsis3 / copernicus - dem - 30 m /
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