On demand mapping services

The most familiar mapping services are those we often use to navigate with our phones for example: google maps. The difference between google maps and many of the examples we have looked at so far is that we do not download data in advance, process it and then create a map. Instead we have access to data when we need it and a a relevant resolution. If we want to make a global map that is 1000 pixels across then Australia is probably just a few dozen pixels wide. If we instead decide to make a map to navigate across Melbourne, then a much higher resolution is needed but only within the boundary of the map.

Cartopy provides access to various of the online mapping services that will serve image data on demand in the form of small image tiles at a specified resolution. The tools automatically query the service and assemble the tiles to make the map but there are some tricks that we need to know before we can use them.

%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.img_tiles as cimgt


# from matplotlib.transforms import offset_copy

from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker

Let’s pre-define some features that we can add to maps later. This can be especially helpful for satellite images where it can be hard to orient the images without some overlays.

coastline = cfeature.NaturalEarthFeature('physical', 'coastline', '10m',
                           edgecolor=(0.0,0.0,0.0),
                           facecolor="none")

rivers = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '10m',
                                        edgecolor='Blue', facecolor="none")


coastline50 = cfeature.NaturalEarthFeature('physical', 'coastline', '50m',
                           edgecolor=(1.0,1.0,1.0),
                           facecolor="none")

rivers50 = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '50m',
                                        edgecolor='Blue', facecolor="none")


lakes = cfeature.NaturalEarthFeature('physical', 'lakes', '10m',
                                        edgecolor="blue", facecolor="blue")

ocean = cfeature.NaturalEarthFeature('physical', 'ocean', '10m',
                           edgecolor="green",
                           facecolor="blue")

## The full list of available interfaces is found in the source code for this one:
## https://github.com/SciTools/cartopy/blob/master/lib/cartopy/io/img_tiles.py

## Continental US terrain images (these no longer work without an api key)
# stamen_Terrain = cimgt.Stamen('terrain-background')
# stamen_TerrainPlus = cimgt.Stamen('terrain')
# stamen_Artist = cimgt.Stamen('watercolor')


## Mapquest satellite / streetmap images 
map_quest_aerial = cimgt.MapQuestOpenAerial()
map_quest_street = cimgt.MapQuestOSM()

## Open Street map
open_street_map = cimgt.OSM()

## Mapbox Satellite images 

mapbox_satellite = cimgt.MapboxTiles(map_id='satellite', 
                                     access_token='pk.eyJ1IjoibG91aXNtb3Jlc2kiLCJhIjoiY2pzeG1mZzFqMG5sZDQ0czF5YzY1NmZ4cSJ9.lpsUzmLasydBlS0IOqe5JA')


mapbox_streets = cimgt.MapboxTiles(map_id='streets', 
                                     access_token='pk.eyJ1IjoibG91aXNtb3Jlc2kiLCJhIjoiY2pzeG1mZzFqMG5sZDQ0czF5YzY1NmZ4cSJ9.lpsUzmLasydBlS0IOqe5JA')


## Google maps image tiles ()
google_maps_street = cimgt.GoogleTiles(style="street") 
google_maps_satellite = cimgt.GoogleTiles(style="satellite") 
google_maps_terrain = cimgt.GoogleTiles(style="terrain") 

Now let’s specify a region of interest using latitude and longitude.

Let’s pick Southern California.

Define the latitude and longitude region (a box) of the area in order to make the map.

# Specify a region of interest

lat0 =  33.5  ; lat1 = 34.5
lon0 =  -119; lon1 = -117.5

map_extent = [lon0, lon1, lat0, lat1]
map_tiles = google_maps_satellite

fig = plt.figure(figsize=(8, 8), facecolor="none")

# Create a GeoAxes in the tile's projection.
ax = plt.axes(projection=map_tiles.crs)
ax.set_extent(map_extent)

# Add the coastline to see where we are 
ax.add_feature(coastline, linewidth=1.5,  edgecolor="Black", zorder=1, alpha=0.5)

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb269c87820>

Here you can set up the size of the map and add on the on-demand image. Try to adjust the resolution and how it works best with the map area.

map_tiles = google_maps_satellite

fig = plt.figure(figsize=(8, 8), facecolor="none")
ax = plt.axes(projection=map_tiles.crs)

# Set the size of the map
ax.set_extent(map_extent)
# Add the on-demand image - the second argument is the resolution and needs to be balanced with the 
# size of the area the map covers. 

ax.add_image(map_tiles, 8, alpha=1.0, zorder=0)
ax.add_feature(coastline, linewidth=1.5,  edgecolor="Black", zorder=1, alpha=0.5)
ax.add_feature(rivers,    linewidth=1.0,  edgecolor="Blue",  zorder=2, alpha=0.5)

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb269a521d0>

/srv/conda/envs/notebook/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_rivers_lake_centerlines.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

Now try to create a map using terrain images rather than satelite imagery

# Specify a region of interest

lat0 =  33.5  ; lat1 = 34.5
lon0 =  -119; lon1 = -117.5

map_extent = [lon0, lon1, lat0, lat1]

map_tiles = google_maps_terrain

fig = plt.figure(figsize=(12, 12), facecolor="none")
# ax = plt.axes(projection=ccrs.PlateCarree(), extent=himalaya_extent)

# Create a GeoAxes in the tile's projection.
ax = plt.axes(projection=map_tiles.crs)

# Set the size of the map
ax.set_extent(map_extent)
# Add the on-demand image - the second argument is the resolution and needs to be balanced with the 
# size of the area the map covers. 

ax.add_image(map_tiles, 8, alpha=1.0, zorder=0)
ax.add_feature(coastline, linewidth=1.5,  edgecolor="Black", zorder=1, alpha=0.5)
ax.add_feature(rivers,    linewidth=1.0,  edgecolor="Blue",  zorder=2, alpha=0.5)

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb2694d5270>

Now save the figure as a PNG file

fig.savefig("LaBasinTerrain.png", dpi=300)

Now let’s try a different region - such as all of Australia. You determine that latitude and longitude bounds to create a map of Australia.

# Specify a region of interest - let's try Australia.

lat0 =    ; lat1 = 
lon0 =    ; lon1 = 

map_extent = [lon0, lon1, lat0, lat1]

map_tiles = google_maps_street

fig = plt.figure(figsize=(8, 8), facecolor="none")

ax = plt.axes(projection=map_tiles.crs)

# Set the size of the map
ax.set_extent(map_extent)

# Add the on-demand image - the second argument is the resolution and needs to be balanced with the 
# size of the area the map covers. 

ax.add_image(map_tiles, 5, alpha=1.0, zorder=2)

ax.add_feature(coastline50, linewidth=1.5,  edgecolor="Black", zorder=1)
ax.add_feature(rivers50,    linewidth=1.0,  edgecolor="Blue",  zorder=2)

  Cell In[15], line 3
    lat0 =    ; lat1 =
              ^
SyntaxError: invalid syntax

## Mapquest - Now let's try mapping the UK

map_extent = [lon0, lon1, lat0, lat1]


fig  = plt.figure(figsize=(8, 8), facecolor="none")
ax1  = plt.subplot(121, projection=google_maps_satellite.crs)
ax2  = plt.subplot(122, projection=google_maps_street.crs)

ax1.set_extent(map_extent)
ax2.set_extent(map_extent)

ax1.add_image(google_maps_satellite, 7, alpha=1.0, zorder=2)
ax2.add_image(google_maps_street, 7, alpha=1.0, zorder=2)

Test the different terrains

map_extent = [lon0, lon1, lat0, lat1]

fig  = plt.figure(figsize=(12, 12), facecolor="none")
ax1  = plt.subplot(121, projection=google_maps_terrain.crs)
ax2  = plt.subplot(122, projection=google_maps_terrain.crs)

ax1.set_extent(map_extent)
ax2.set_extent(map_extent)

ax1.add_image(google_maps_satellite,      7,  alpha=1.0, zorder=2)
ax2.add_image(open_street_map,       7,  alpha=1.0, zorder=2)

# To test what the influence of the resolution parameter is on the map that is returned ... 

lat0 =   51  ; lat1 = 52
lon0 =  -1  ; lon1 = 1

map_extent = [lon0, lon1, lat0, lat1]


##maptype1 = mapbox_satellite
##maptype2 = mapbox_streets

## Could also try these ... 

maptype1 = open_street_map
maptype2 = open_street_map

fig  = plt.figure(figsize=(12, 12), facecolor="none")
ax1  = plt.subplot(221, projection=maptype1.crs)
ax2  = plt.subplot(222, projection=maptype2.crs)
ax3  = plt.subplot(223, projection=maptype1.crs)
ax4  = plt.subplot(224, projection=maptype2.crs)

ax1.set_extent(map_extent)
ax2.set_extent(map_extent)
ax3.set_extent(map_extent)
ax4.set_extent(map_extent)

ax1.add_image(maptype1, 8, alpha=1.0, zorder=2)
ax2.add_image(maptype2, 8, alpha=1.0, zorder=2)

ax3.add_image(maptype1, 10, alpha=1.0, zorder=2)
ax4.add_image(maptype2, 10, alpha=1.0, zorder=2)

fig.savefig("ResolutionTest.png", dpi=600)

See how well that worked by opening the full-size image.

There really is a lot of detail available, but you can also see how long it takes to download and build the different resolutions. There is something of an art to finding the right balance.

You may also find that very high resolution output results in memory errors. If so, restart the kernel and try again with a small image or lower dpi.

## For fun, this is the street where I grew up ... 

lat0 =   49.9 ; lat1 = 49.92
lon0 =   -119.50   ; lon1 = -119.48

map_extent = [lon0, lon1, lat0, lat1]


figE18  = plt.figure(figsize=(12, 6), facecolor="none")
ax1  = plt.subplot(121, projection=google_maps_street.crs)
ax2  = plt.subplot(122, projection=google_maps_satellite.crs)

ax1.set_extent(map_extent)
ax2.set_extent(map_extent)

ax1.add_image(google_maps_street, 17,      alpha=1.0, zorder=2)
ax2.add_image(google_maps_satellite  , 17, alpha=1.0, zorder=2)

figE18.savefig("home.png", dpi=900)

Again, open the full-size image to see how well that worked

Now you try to map the street where you grew up!

# Try it here!