目录

上一个主题

11.1. Cartopy 介绍

下一个主题

11.3. Cartopy绘图要素

关注公众号


常见问题

  1. Windows下的安装说明
  2. Jupyter免费在线实验环境
  3. 勘误与补充


11.2. 开始使用 Cartopy

11.2.1. 自定义边界形状

此示例演示如何使用自定义形状几何图形而不是投影的默认边界。 在本例中,我们将边界定义为轴坐标中的圆。这意味着无论地图本身的范围如何,边界始终是一个圆。

>>> %matplotlib inline
>>>
>>> import matplotlib.path as mpath
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>>
>>> import cartopy.crs as ccrs
>>> import cartopy.feature as cfeature
>>>
>>>
>>> fig = plt.figure(figsize=[10, 5])
>>> ax1 = fig.add_subplot(1, 2, 1, projection=ccrs.SouthPolarStereo())
>>> ax2 = fig.add_subplot(1, 2, 2, projection=ccrs.SouthPolarStereo(),
>>>                       sharex=ax1, sharey=ax1)
>>> fig.subplots_adjust(bottom=0.05, top=0.95,
>>>                     left=0.04, right=0.95, wspace=0.02)
>>>
>>> # Limit the map to -60 degrees latitude and below.
>>> ax1.set_extent([-180, 180, -90, -60], ccrs.PlateCarree())
>>>
>>> ax1.add_feature(cfeature.LAND)
>>> ax1.add_feature(cfeature.OCEAN)
>>>
>>> ax1.gridlines()
>>> ax2.gridlines()
>>>
>>> ax2.add_feature(cfeature.LAND)
>>> ax2.add_feature(cfeature.OCEAN)
>>>
>>> # Compute a circle in axes coordinates, which we can use as a boundary
>>> # for the map. We can pan/zoom as much as we like - the boundary will be
>>> # permanently circular.
>>> theta = np.linspace(0, 2*np.pi, 100)
>>> center, radius = [0.5, 0.5], 0.5
>>> verts = np.vstack([np.sin(theta), np.cos(theta)]).T
>>> circle = mpath.Path(verts * radius + center)
>>>
>>> ax2.set_boundary(circle, transform=ax2.transAxes)
>>>
>>> plt.show()
>>>
_images/cartopy-begin_1_0.png

11.2.2. 卡特里娜飓风

这个例子使用Shapely的功能来说明可能受到卡特里娜飓风严重影响的州。

>>> import matplotlib.patches as mpatches
>>> import matplotlib.pyplot as plt
>>> import shapely.geometry as sgeom
>>>
>>> import cartopy.crs as ccrs
>>> import cartopy.io.shapereader as shpreader
>>>
>>>
>>> def sample_data():
>>>     """
>>>     Return a list of latitudes and a list of longitudes (lons, lats)
>>>     for Hurricane Katrina (2005).
>>>
>>>     The data was originally sourced from the HURDAT2 dataset from AOML/NOAA:
>>>     http://www.aoml.noaa.gov/hrd/hurdat/newhurdat-all.html on 14th Dec 2012.
>>>
>>>     """
>>>     lons = [-75.1, -75.7, -76.2, -76.5, -76.9, -77.7, -78.4, -79.0,
>>>             -79.6, -80.1, -80.3, -81.3, -82.0, -82.6, -83.3, -84.0,
>>>             -84.7, -85.3, -85.9, -86.7, -87.7, -88.6, -89.2, -89.6,
>>>             -89.6, -89.6, -89.6, -89.6, -89.1, -88.6, -88.0, -87.0,
>>>             -85.3, -82.9]
>>>
>>>     lats = [23.1, 23.4, 23.8, 24.5, 25.4, 26.0, 26.1, 26.2, 26.2, 26.0,
>>>             25.9, 25.4, 25.1, 24.9, 24.6, 24.4, 24.4, 24.5, 24.8, 25.2,
>>>             25.7, 26.3, 27.2, 28.2, 29.3, 29.5, 30.2, 31.1, 32.6, 34.1,
>>>             35.6, 37.0, 38.6, 40.1]
>>>
>>>     return lons, lats
>>>
>>>
>>> fig = plt.figure()
>>> ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.LambertConformal())
>>>
>>> ax.set_extent([-125, -66.5, 20, 50], ccrs.Geodetic())
>>>
>>> shapename = 'admin_1_states_provinces_lakes_shp'
>>> states_shp = shpreader.natural_earth(resolution='110m',
>>>                                      category='cultural', name=shapename)
>>>
>>> lons, lats = sample_data()
>>>
>>> # to get the effect of having just the states without a map "background"
>>> # turn off the outline and background patches
>>> ax.background_patch.set_visible(False)
>>> ax.outline_patch.set_visible(False)
>>>
>>> ax.set_title('US States which intersect the track of '
>>>              'Hurricane Katrina (2005)')
>>>
>>> # turn the lons and lats into a shapely LineString
>>> track = sgeom.LineString(zip(lons, lats))
>>>
>>> # buffer the linestring by two degrees (note: this is a non-physical
>>> # distance)
>>> track_buffer = track.buffer(2)
>>>
>>> def colorize_state(geometry):
>>>     facecolor = (0.9375, 0.9375, 0.859375)
>>>     if geometry.intersects(track):
>>>         facecolor = 'red'
>>>     elif geometry.intersects(track_buffer):
>>>         facecolor = '#FF7E00'
>>>     return {'facecolor': facecolor, 'edgecolor': 'black'}
>>>
>>> ax.add_geometries(
>>>     shpreader.Reader(states_shp).geometries(),
>>>     ccrs.PlateCarree(),
>>>     styler=colorize_state)
>>>
>>> ax.add_geometries([track_buffer], ccrs.PlateCarree(),
>>>                   facecolor='#C8A2C8', alpha=0.5)
>>> ax.add_geometries([track], ccrs.PlateCarree(),
>>>                   facecolor='none', edgecolor='k')
>>>
>>> # make two proxy artists to add to a legend
>>> direct_hit = mpatches.Rectangle((0, 0), 1, 1, facecolor="red")
>>> within_2_deg = mpatches.Rectangle((0, 0), 1, 1, facecolor="#FF7E00")
>>> labels = ['State directly intersects\nwith track',
>>>           'State is within \n2 degrees of track']
>>> ax.legend([direct_hit, within_2_deg], labels,
>>>           loc='lower left', bbox_to_anchor=(0.025, -0.1), fancybox=True)
>>>
>>> plt.show()
>>>
/usr/lib/python3/dist-packages/cartopy/io/__init__.py:260: DownloadWarning: Downloading: http://naciscdn.org/naturalearth/110m/cultural/ne_110m_admin_1_states_provinces_lakes_shp.zip
  warnings.warn('Downloading: {}'.format(url), DownloadWarning)
_images/cartopy-begin_3_1.png

11.2.3. 使用Cartopy和AxesGrid工具包

此示例演示如何将 Cartopy GeoAxes 与来自 mpl_toolkits.axes_grid1 的 AxesGrid 一起使用。 脚本使用Plate Carree投影构造axes_类kwarg并将其传递给AxesGrid实例。 AxesGrid内置标签关闭,而是使用创建网格线的标准过程。然后绘制一些假数据。

>>> import cartopy.crs as ccrs
>>> from cartopy.mpl.geoaxes import GeoAxes
>>> from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.axes_grid1 import AxesGrid
>>> import numpy as np
>>>
>>>
>>> def sample_data_3d(shape):
>>>     """Return `lons`, `lats`, `times` and fake `data`"""
>>>     ntimes, nlats, nlons = shape
>>>     lats = np.linspace(-np.pi / 2, np.pi / 2, nlats)
>>>     lons = np.linspace(0, 2 * np.pi, nlons)
>>>     lons, lats = np.meshgrid(lons, lats)
>>>     wave = 0.75 * (np.sin(2 * lats) ** 8) * np.cos(4 * lons)
>>>     mean = 0.5 * np.cos(2 * lats) * ((np.sin(2 * lats)) ** 2 + 2)
>>>
>>>     lats = np.rad2deg(lats)
>>>     lons = np.rad2deg(lons)
>>>     data = wave + mean
>>>
>>>     times = np.linspace(-1, 1, ntimes)
>>>     new_shape = data.shape + (ntimes, )
>>>     data = np.rollaxis(data.repeat(ntimes).reshape(new_shape), -1)
>>>     data *= times[:, np.newaxis, np.newaxis]
>>>
>>>     return lons, lats, times, data
>>>
>>>
>>>
>>> projection = ccrs.PlateCarree()
>>> axes_class = (GeoAxes,
>>>               dict(map_projection=projection))
>>>
>>> lons, lats, times, data = sample_data_3d((6, 73, 145))
>>>
>>> fig = plt.figure()
>>> axgr = AxesGrid(fig, 111, axes_class=axes_class,
>>>                 nrows_ncols=(3, 2),
>>>                 axes_pad=0.6,
>>>                 cbar_location='right',
>>>                 cbar_mode='single',
>>>                 cbar_pad=0.2,
>>>                 cbar_size='3%',
>>>                 label_mode='')  # note the empty label_mode
>>>
>>> for i, ax in enumerate(axgr):
>>>     ax.coastlines()
>>>     ax.set_xticks(np.linspace(-180, 180, 5), crs=projection)
>>>     ax.set_yticks(np.linspace(-90, 90, 5), crs=projection)
>>>     lon_formatter = LongitudeFormatter(zero_direction_label=True)
>>>     lat_formatter = LatitudeFormatter()
>>>     ax.xaxis.set_major_formatter(lon_formatter)
>>>     ax.yaxis.set_major_formatter(lat_formatter)
>>>
>>>     p = ax.contourf(lons, lats, data[i, ...],
>>>                     transform=projection,
>>>                     cmap='RdBu')
>>>
>>> axgr.cbar_axes[0].colorbar(p)
>>>
>>> plt.show()
>>>
_images/cartopy-begin_5_0.png

11.2.4. Cartopy标志

生成 Cartopy 徽标的实际代码。

>>> import cartopy.crs as ccrs
>>> import matplotlib.pyplot as plt
>>> import matplotlib.textpath
>>> import matplotlib.patches
>>> from matplotlib.font_manager import FontProperties
>>> import numpy as np
>>>
>>> fig = plt.figure(figsize=[12, 6])
>>> ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson())
>>>
>>> ax.coastlines()
>>> ax.gridlines()
>>>
>>> # generate a matplotlib path representing the word "cartopy"
>>> fp = FontProperties(family='Bitstream Vera Sans', weight='bold')
>>> logo_path = matplotlib.textpath.TextPath((-175, -35), 'cartopy',
>>>                                          size=1, prop=fp)
>>> # scale the letters up to sensible longitude and latitude sizes
>>> logo_path._vertices *= np.array([80, 160])
>>>
>>> # add a background image
>>> im = ax.stock_img()
>>> # clip the image according to the logo_path. mpl v1.2.0 does not support
>>> # the transform API that cartopy makes use of, so we have to convert the
>>> # projection into a transform manually
>>> plate_carree_transform = ccrs.PlateCarree()._as_mpl_transform(ax)
>>> im.set_clip_path(logo_path, transform=plate_carree_transform)
>>>
>>> # add the path as a patch, drawing black outlines around the text
>>> patch = matplotlib.patches.PathPatch(logo_path,
>>>                                      facecolor='none', edgecolor='black',
>>>                                      transform=ccrs.PlateCarree())
>>> ax.add_patch(patch)
>>>
>>> plt.show()
_images/cartopy-begin_7_0.png