Maps

1. Introduction

Maps are a way to present information on a (roughly) spherical earth on a flat plane, like a page or a screen. Here are two examples of common map projections. The projection is only accurate in the region where the plane touches the sphere, and is less accurate as the distance between the plane and the sphere increases.

Mercator

Lambert conformal conic

You can read more about map projections from Map Projections – a Working Manual, the source of the images above, or, more entertainingly, from XKCD.

We'll use cartopy to plot on maps. Check out the gallery for inspiration.

Here we have the most basic projection: plate carrée, which is an equirectangular projection, and is essentially equivalent to just plotting the longitude and latitude values without a projection. I will refer to longitude and latitude as "geographic coordinates".

We can make an axes that is plotted in geographic coordinates (or, indeed, any projection we choose) by using the projection keyword argument to fig.add_subplot(). Here we also plot the coastline and add gridlines.

cartopy provides a number of projections. Available projections are:

PlateCarree
AlbersEqualArea
AzimuthalEquidistant
LambertConformal
LambertCylindrical
Mercator
Miller
Mollweide
Orthographic
Robinson
Sinusoidal
Stereographic
TransverseMercator
UTM
InterruptedGoodeHomolosine
RotatedPole
OSGB
EuroPP
Geostationary
Gnomonic
LambertAzimuthalEqualArea
NorthPolarStereo
OSNI
SouthPolarStereo  

Lambert Conformal Conic is a useful projection in numerical modeling because it preserves right angles. Here we use the projection without any keyword specifications, but with the coastline plotted so that we have something to look at.

The projection that we choose in the axes line with projection= is the projection that the plot is in. Data from any projection can be plotted on this map, but we will have to tell it which projection it is in.

Let's make a map of the Gulf of Mexico using the LambertConformal projection. Projections take in different keywords to specify properties. For this projection, we can specify the central longitude and latitude, which control the center of the projection. Our selection in the example is not far off from the default, so it looks similar to the previous plot.

Don't forget to center your projections around your area of interest to be as accurate as possible for your purposes.

The map is plotted in a projected coordinate system, with units in meters, but the package deals with the projection behind the scenes. We can see this by looking at the limits of the two axes, which don't look like longitude/latitude at all:

This same call to a plot set up with the PlateCarree projection, which is in geographic coordinates (lon/lat) does give limits in longitude and latitude, because in that case we told the plot to be in those coordinates, but in this case we said to use a LambertConformal.

We can use whatever type of coordinates we want, including latitude and longitude, as long as we tell cartopy which type we are using.

As you saw above, we set the limits of the plot not with xlim and ylim, but with extent and the appropriate projection object.

Exercise

Create a map of the Gulf of Mexico using a different projection. How does it compare to the map above?

This is pretty good, but there are some limitations in this package currently. One is that we can't add labels to the lat/lon lines for the Lambert Conformal Conic projection. We can do this using Mercator, though:

When we want to add something to the plot, we just need to tell it what projection the information is given in using the transform keyword argument.

If the information is in latitude/longitude – typical for the way people tend to think about information (instead of projected locations) – then we give the Plate Carree projection with the transform keyword argument to the plot call:

transform=ccrs.PlateCarree()

For example, to plot some points with a particular projection, you can type:

plt.plot(xpts, ypts, transform=ccrs.projection_that_xpts_and_ypts_are_given_in)

A nice thing about the cartopy package is that you can plot directly data from any projection — you just tell it the projection through the transform keyword argument when you add to the plot.

Exercise

The latitude and longitude of College Station are given below. Plot the location of College Station on the map above with a red dot.

lat_cll = 30.0 + 36.0/60.0 + 5.0/3600.0
lon_cll = -(96.0 + 18.0/60.0 + 52.0/3600.0)

What happens if you put in the wrong projection or no projection?

Exercise

Data from any projection can be added to a map, the data must just be input with its projection using the transform keyword.

The x, y location of Austin, TX, is given below in the Mercator projection. Plot the location of Austin in Mercator coordinates on the map above with a blue 'x'.

x, y = -10880707.173023093, 3516376.324225941

Point conversion

While cartopy removes the need to convert points on your own between projections (instead doing it behind the scenes), you can always convert between projections if you want using the following. Or, if you want to transform more than one point, use projection.transform_points(projection, x, y).

Exercise

Convert the Mercator coordinates given for Austin to latitude and longitude and confirm that they are correct.

Other features you can add

The code we used earlier, like:

ax.add_feature(cartopy.feature.LAND)

was a convenience function wrapping more complex and capable code options. Here we explore a little more the capabilities. Note that this requires downloading data which you will see a warning about the first time you run the code.

We can set up the ability to plot with high resolution land data:

There are also some built-in colors, but you can use any matplotlib color available to color the land or water.

Using higher resolution can be pretty significantly different.

Here we will prepare the higher resolution land and ocean information for the highest resolution available, then use it in the plot.

Here is a list (with reference names in some cases appended) of the many features that are available through Natural Earth:

(10, 50, 110 for high, medium, low resolution)

Physical Vector Data Themes:

(physical)

• Coastline (10, 50, 110): coastline
• Land (10, 50, 110): land
• Ocean (10, 50, 110): ocean
• Minor Islands (10): minor_islands, minor_islands_coastline
• Reefs (10): reefs
• Physical region features (10): geography_regions_polys, geography_regions_points, geography_regions_elevation_points, geography_marine_polys
• Rivers and Lake Centerlines (10, 50, 110): rivers_lake_centerlines
• Lakes (10, 50, 110): lakes
• Glaciated areas (10, 50, 110): glaciated_areas
• Antarctic ice shelves (10, 50): antarctic_ice_shelves_polys, antarctic_ice_shelves_lines
• Bathymetry (10): bathymetry_all or choose which depth(s)
• Geographic lines (10, 50): geographic_lines
• Graticules (10, 50, 110): (grid lines) graticules_all or choose degree interval

Raster Data Themes:

(raster: land coloring)

• Cross Blended Hypsometric Tints (10, 50)
• Natural Earth 1 (10, 50)
• Natural Earth 2 (10, 50)
• Ocean Bottom (10, 50)
• Bathymetry (50)
• Shaded Relief (10, 50)
• Gray Earth (10, 50)
• Manual Shaded Relief (10, 50)

Cultural Vector Data Themes:

(cultural)

• Countries (10, 50, 110): admin_0_countries, admin_0_countries_lakes, admin_0_boundary_lines
• Disputed areas and breakaway regions (10, 50)
• First order admin (provinces, departments, states, etc.) (10, 50): e.g. admin_1_states_provinces_lines
• Populated places (10, 50, 110)
• Urban polygons (10, 50)
• Parks and protected areas (10): parks_and_protected_lands
• Pacific nation groupings (10, 50, 110)
• Water boundary indicators (10)

Here is an example showing state boundaries: