Plotting in Cartopy

I have been using Basemap for years in my Python scripting. Now I would like to slowly switch to Python 3, but Basemap is no longer officially supported there. There seems to be some workaround (such as using customized versions of Basemap), but personally I do not want to try that. So, I guess it is time to find an alternative, such as Cartopy.

In this post, I leave a record of my learning more about Cartopy, as well as some useful tricks (that I found out after hours of experiments..)

1. a basic map with raster data

import numpy as np
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
fig1 = plt.figure(figsize=(9,7))

# setting up projection
ax1 = plt.subplot(1,1,1, projection=ccrs.Mercator())

# add basic layers
ax1.add_feature(cartopy.feature.OCEAN, linewidth=1, facecolor='lightgrey', edgecolor='k', zorder=0)
ax1.add_feature(cartopy.feature.LAND, linewidth=1, facecolor='cornsilk', edgecolor='k', zorder=1)
ax1.add_feature(cartopy.feature.LAKES, linewidth=1, facecolor='lightgrey', edgecolor='k', zorder=2)

# add state boundaries
state = cartopy.feature.NaturalEarthFeature(category='cultural', scale='50m', facecolor='none', \
                                            name='admin_1_states_provinces_shp')
ax1.add_feature(state, zorder=11)

# add data
ax1.pcolor(lons, lats, np.ma.masked_array(ERA_map, mask=CONUS_mask), \
           cmap=cmap_3color, norm=norm_3color, transform=ccrs.PlateCarree(), zorder=10)

# set map extent
ax1.set_extent([-130, -63, 23, 46], crs=ccrs.Geodetic())
plt.show()

Here is my take to get Cartopy to work:

 

  1. Projection needs to be defined when creating the plot (or subplot).
  2. Ocean, land, lakes are added through “ax.add_feature” function. Data is from cartopy.feature.XXXX. When called for the first time, it might take some time to download the data, so just be patient.
  3. Country boundaries are added through “ax.add_feature” function. Data is in the following format: cartopy.feature.NaturalEarthFeature(category=<arg1>, scale=<arg2>, facecolor=<arg3>, name=<arg4>)

Also a note on the data transforamtion as in “pcolor” command. With some testing, it turns out only ccrs.PlateCarree() projection would work, but it works with maps of any projection (as defined in the axes creation). So this can be fixed.

 

2. Switching between projections

Since I am frequently working with WRF, I need to handle data in Lambert Conformal projection. Switching between this Lambert projection and lat/lon coordinate (i.e. Geodetic) turns out to be quite simple. Just take this command:

import cartopy.crs as ccrs
lccproj = ccrs.LambertConformal(central_longitude=ref_lon, central_latitude=ref_lat, 
                                standard_parallels=(par_lat1, par_lat2), globe=None, cutoff=10)
lcc_x, lcc_y = lccproj.transform_point(lon, lat, latlonproj)
lon, lat = latlonproj.transform_point(lcc_x, lcc_y, lccproj)

lcc_x and lcc_y are the coordinates in Lambert projection, while lat and lon are the latitude and longitude of this point.

Just a comparison, in Basemap we use:

from mpl_toolkits.basemap import Basemap
m = Basemap(projection='lcc', .......)
lcc_x, lcc_y = m(lon, lat)
lon, lat = m(lcc_x, lcc_y, inverse=True)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s