Georeference a Static Image API

index.html

"""static.py

An example of fetching an image from the Mapbox Static Image API with a
matching "world file" for georeferencing the image.
"""

import os

import mercantile
import requests

# Parameters describing the static map.
center_lon = -109.124
center_lat = 36.9736
zoom_lvl = 6
width = 1024
height = 1024
access_token = os.getenv("MapboxAccessToken")

# Construct a URL for the static map.
url = "https://api.mapbox.com/styles/v1/mapbox/satellite-v9/static/{:f},{:f},{:d},0/{:d}x{:d}".format(
    center_lon, center_lat, zoom_lvl, width, height
)

# Calculate the georeferencing parameters for the static map.  The
# coordinate reference system for static map images is spherical web
# mercator, known to GIS systems as "EPSG:3857".  Methods of the Python
# mercantile package are needed to go back and forth between EPSG:3857
# and geographic longitude and latitude.
#
# First, we compute the resolution (in meters per pixel) of the static
# map image.
center_tile = mercantile.tile(center_lon, center_lat, zoom_lvl)
left, bottom, right, top = mercantile.xy_bounds(center_tile)
xres = (right - left) / 512
yres = (bottom - top) / 512

# Next we compute the EPSG:3857 coordinates of the center of the upper
# leftmost pixel in the static map image.
center_x, center_y = mercantile.xy(center_lon, center_lat)
image_ul_cx = center_x - (xres / 2) * (width - 1)
image_ul_cy = center_y - (yres / 2) * (height - 1)

# Finally, create the content of the world file. It's a text file with 6
# lines. The first line is the change in coordinate value per pixel in
# the x direction (left to right in the image). The second and third
# lines parametrize rotation and are usually 0.0. The fourth line in the
# change in coordinate value per pixel in the y direction (top to
# bottom). The fifth and sixth lines are the x and y coordinates of the
# center of the image's upper leftmost pixel.
world_file_content = """{:f}
0
0
{:f}
{:f}
{:f}
""".format(
    xres, yres, image_ul_cx, image_ul_cy
)

# Then write out the files using the same basename and these required
# file extensions.
with open("static.jpw", "w") as fjpw:
    fjpw.write(world_file_content)

with open("static.jpg", "wb") as fjpg:
    resp = requests.get(url, params=dict(access_token=access_token))
    fjpg.write(resp.content)



Let's say we've extracted a wedge shape from the static map image and the shape has corners with image x,y coordinates (with origin at the upper left) like

feature_coords = [[10, 20], [100, 45], [80, 55], [10, 20]]

Our Python rasterio package has a handy API for transforming these image coordinates back to world coordinates.

with rasterio.open("static.jpg") as dataset:
    geotransform = dataset.transform

mercator_coords = [geotransform * (x, y) for (x, y) in feature_coords]
# mercator_coords: [(-12761570.324507501, 5037139.8929225), (-12651501.0037375, 5006565.0815975005), (-12675960.8527975, 4994335.1570675), (-12761570.324507501, 5037139.8929225)]
And then a method of the mercantile package can be used to go back to lon/lat.

lonlat_coords = [mercantile.lnglat(x, y) for (x, y) in mercator_coords]
# lonlat_coords: [LngLat(lng=-114.6391367187121, lat=41.16790936699352), LngLat(lng=-113.65036718710863, lat=40.96082503555483), LngLat(lng=-113.87009374968719, lat=40.87780876755107), LngLat(lng=-114.6391367187121, lat=41.16790936699352)]
(Edited to fix the coordinate values)