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:
- Projection needs to be defined when creating the plot (or subplot).
- 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.
- 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)