Joseph Barraud
This notebook explains how to make maps in orthographic projection with matplotlib basemap
and cartopy
. The interpies
library is also used to load the data and create images. interpies
is available here.
Note: For this notebook to work, you need to install interpies
and its dependencies, as well as basemap
and cartopy
. Using conda
, that should be as easy as writing:
conda install basemap
conda install cartopy
The orthographic projection is typically the one that results in an image of the Earth as a globe floating in space. Provided you have the right image, it is relatively easy to warp it around the globe and make a map that looks both interesting and unusual.
Let's start by importing the interpies
module, basemap
, cartopy
and a few others.
import sys, os
import numpy as np
import imageio
from matplotlib import pyplot as plt
# basemap
from mpl_toolkits.basemap import Basemap
# cartopy
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
# interpies
import interpies
interpies.__version__
'0.3.1'
# plot within the notebook
%matplotlib inline
fig, ax = plt.subplots(figsize=(10,10))
m = Basemap(projection='ortho', lat_0=45, lon_0=0, resolution='c', ax=ax)
# draw coastlines and borders
m.drawcoastlines()
m.drawcountries()
# draw meridians and parallels
m.drawmeridians(np.arange(0,360,30))
m.drawparallels(np.arange(-90,90,30))
plt.show()
C:\Anaconda3\envs\interpies_test\lib\site-packages\mpl_toolkits\basemap\__init__.py:1711: MatplotlibDeprecationWarning: The axesPatch function was deprecated in version 2.1. Use Axes.patch instead. if limb is not ax.axesPatch:
cartopy
¶Similarly, the orientation of the globe is controlled by parameters of the projection class.
Regarding meridians and parallels, the defaults are different but I tried to reproduce the behaviour of basemap
here.
fig = plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic(0, 45))
# draw coastlines and borders
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, lw=0.5)
# draw meridians and parallels
gl = ax.gridlines(color='k', linestyle=(0, (1, 1)), xlocs=range(0,390,30), ylocs=[-80, -60, -30, 0, 30, 60, 80])
fig, ax = plt.subplots(figsize=(10, 10))
m = Basemap(projection='ortho', lat_0=45, lon_0=0, resolution='l', ax=ax)
m.etopo()
# draw borders only
m.drawcountries()
m.drawmeridians(np.arange(0,360,30))
m.drawparallels(np.arange(-90,90,30))
plt.show()
C:\Anaconda3\envs\interpies_test\lib\site-packages\mpl_toolkits\basemap\__init__.py:1711: MatplotlibDeprecationWarning: The axesPatch function was deprecated in version 2.1. Use Axes.patch instead. if limb is not ax.axesPatch:
cartopy
makes use of Natural Earth data for its default background image.
fig = plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic(0, 45))
# add shaded relief image
ax.stock_img()
# draw coastlines and borders
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, lw=0.5)
# draw meridians and parallels
gl = ax.gridlines(color='k', linestyle=(0, (1, 1)), xlocs=range(0,390,30), ylocs=[-80, -60, -30, 0, 30, 60, 80])
Now we want to display our own data in this orthographic projection. We need a global dataset, so let's try the satellite gravity anomaly map from the Scripps Institution of Oceanography. The dataset is commonly known as Sandwell's marine gravity data (Sandwell et al., 2014, also NOAA and NGA).
In basemap
, the function behind the etopo
and bluemarble
methods is warpimage
: it takes an image and distorts it to fit the geometry of the map. In cartopy
, the equivalent is the imshow
method.
In both cases, the first thing to do is to create an image of the gravity anomaly map. As the original data is provided as a map in a Mercator projection, it needs to be "unprojected". The method I used to create a 5x5 min grid can be found here. However, since we don't need the full resolution for a global map, I resampled the grid further to a 10 minute cell size.
Using interpies
, first load the data in a Grid
object.
grid1 = interpies.open('../data/sandwell_grav_v24_10min.tif')
Next display the grid with the show
method.
ax = grid1.show(figsize=(18, 12))
With info
, we can see that the grid does not cover all the latitudes from -90 to +90, so we will have to deal with that later.
grid1.info()
* Info * Grid name: sandwell_grav_v24_10min Filename: ../data/sandwell_grav_v24_10min.tif Coordinate reference system: epsg:4326 Grid size: 2160 columns x 969 rows Cell size: 0.1667 x 0.1667 Lower left corner (pixel centre): (0.083, -80.679) Grid extent (outer limits): west: 0.000, east: 360.000, south: -80.762, north: 80.738 No Data Value: -99999.0 Number of null cells: 0 (0.00%) * Statistics * mean = -0.2993307946330301 sigma = 31.987015187122367 min = -352.09283447265625 max = 699.9725952148438
But first, let's export an image of the map in a png format. Using the save_image
method, one can create an image of the grid with hillshading and colours but without the labels and the colorbar. The function has the same parameters as show
. Here, I am using cmap_brightness
and hs_contrast
to make the map a little bit brighter.
grid1.save_image('../data/grav_v24_10min.png', cmap_brightness=1.5, hs_contrast=3.0)
The grid was successfully saved as an image in ../data/grav_v24_10min.png
<matplotlib.figure.Figure at 0x240013b7400>
In cartopy
, the imshow
method is used to display images on axes. The input is the RGB image created previously, so it needs to be loaded in an array.
im = imageio.imread('../data/grav_v24_10min.png')
im.shape
(969, 2160, 4)
The image has a fourth channel (alpha) that for some reason will create display problems with imshow
. So it needs to be removed.
The other important parameter that needs tweaking is the extent. Reading from the grid info, it is [west, east, south, north]
. However, imshow
would not accept both west=0
if east=360
, so I added a tiny shift and started the grid at west=0.001
.
fig = plt.figure(figsize=(12, 12))
ax = plt.axes(projection=ccrs.Orthographic(0, 45))
ax.imshow(im[:,:,:3], extent=[0.001, 360, -80.720, 80.738], origin='upper', transform=ccrs.PlateCarree())
# draw coastlines
ax.add_feature(cfeature.COASTLINE)
# draw meridians and parallels
gl = ax.gridlines(color='k', linestyle=(0, (1, 1)), xlocs=range(0,390,30), ylocs=[-80, -60, -30, 0, 30, 60, 80])
Here is another view:
fig = plt.figure(figsize=(12, 12))
ax = plt.axes(projection=ccrs.Orthographic(-60, -15))
ax.imshow(im[:,:,:3], extent=[0.001, 360, -80.720, 80.738], origin='upper', transform=ccrs.PlateCarree())
# draw coastlines
ax.add_feature(cfeature.COASTLINE)
# draw meridians and parallels
gl = ax.gridlines(color='k', linestyle=(0, (1, 1)), xlocs=range(0,390,30), ylocs=[-80, -60, -30, 0, 30, 60, 80])
Things are a little more complicated with basemap
.
From the docstring of warpimage
, we learn that the image needs to cover the entire globe and start at -180W and the South Pole.
In order to make sure the image covers the whole globe, we need to add some padding to fill the gap between the top and bottom edges of the image and the poles. With a cell size of 10 minutes (or 0.16666 degrees), a grid going from -90 to +90 would normally have 180/(10/60)=1080 rows. Our image has 969 rows, so 111 rows are missing. As this is an odd number, we cannot add the same number of rows on both sides. At that resolution, it shouldn't be a problem.
# add 56 rows at the top
im_pad = np.insert(im, 0, 255*np.ones((56, 2160, 4)), axis=0)
# add 55 rows at the bottom
im_pad = np.insert(im_pad, 969+56, 255*np.ones((55, 2160, 4)), axis=0)
Let's see the result:
fig, ax = plt.subplots(figsize=(18,12))
ax.imshow(im_pad);
According to the documentation, the global map should be centred on the Greenwich (0 degree) meridian. We can easily achieve that by using the roll
function:
# roll
im_pad_roll = np.roll(im_pad, shift=(0, 1080, 0), axis=(0, 1, 2))
We can now use warpimage
to display our padded image on the globe. One last thing however: the function requires the image to be read from a file. So let's save the previous image to a png file.
# Save the image
imageio.imwrite('../data/grav_v24_10min_pad_roll.png', im_pad_roll)
fig,ax = plt.subplots(figsize=(12, 12))
m = Basemap(projection='ortho', lat_0=45, lon_0=0, resolution='c', ax=ax)
m.warpimage(image='../data/grav_v24_10min_pad_roll.png')
# draw coastlines.
m.drawcoastlines()
# draw meridians and parallels
m.drawmeridians(np.arange(0, 360, 30))
m.drawparallels(np.arange(-90, 90, 30))
plt.show()
C:\Anaconda3\envs\interpies_test\lib\site-packages\mpl_toolkits\basemap\__init__.py:1711: MatplotlibDeprecationWarning: The axesPatch function was deprecated in version 2.1. Use Axes.patch instead. if limb is not ax.axesPatch:
Sandwell, D. T., R. D. Müller, W. H. F. Smith, E. Garcia, R. Francis. 2014. New global marine gravity model from CryoSat-2 and Jason-1 reveals buried tectonic structure, Science, Vol. 346, no. 6205, pp. 65-67, doi: 10.1126/science.1258213.