# 11.5. Cartopy地图绘图2¶

## 11.5.1. 使用错误空间参考的问题¶

>>> %matplotlib inline
>>>
>>> import cartopy.crs as ccrs
>>> import cartopy.feature as cfeature
>>> from cartopy.io.img_tiles import Stamen
>>> import matplotlib.pyplot as plt
>>> from matplotlib.lines import Line2D as Line
>>> from matplotlib.patheffects import Stroke
>>> import numpy as np
>>> import shapely.geometry as sgeom
>>> from shapely.ops import transform as geom_transform
>>>
>>>
>>> def transform_fn_factory(target_crs, source_crs):
>>>     """
>>>     Return a function which can be used by shapely.op.transform
>>>     to transform the coordinate points of a geometry.
>>>
>>>     The function explicitly *does not* do any interpolation or clever
>>>     transformation of the coordinate points, so there is no guarantee
>>>     that the resulting geometry would make any sense.
>>>
>>>     """
>>>     def transform_fn(x, y, z=None):
>>>         new_coords = target_crs.transform_points(source_crs,
>>>                                                  np.asanyarray(x),
>>>                                                  np.asanyarray(y))
>>>         return new_coords[:, 0], new_coords[:, 1], new_coords[:, 2]
>>>
>>>     return transform_fn
>>>
>>>
>>>
>>> # Define the two coordinate systems with different ellipses.
>>> wgs84 = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84',
>>>                                           ellipse='WGS84'))
>>> sphere = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84',
>>>                                            ellipse='sphere'))
>>>
>>> # Define the coordinate system of the data we have from Natural Earth and
>>> # acquire the 1:10m physical coastline shapefile.
>>> geodetic = ccrs.Geodetic(globe=ccrs.Globe(datum='WGS84'))
>>> dataset = cfeature.NaturalEarthFeature(category='physical',
>>>                                        name='coastline',
>>>                                        scale='10m')
>>>
>>> # Create a Stamen map tiler instance, and use its CRS for the GeoAxes.
>>> tiler = Stamen('terrain-background')
>>> fig = plt.figure()
>>> ax = fig.add_subplot(1, 1, 1, projection=tiler.crs)
>>> ax.set_title('The effect of incorrectly referencing the Solomon Islands')
>>>
>>> # Pick the area of interest. In our case, roughly the Solomon Islands, and
>>> # get hold of the coastlines for that area.
>>> extent = [155, 163, -11.5, -6]
>>> ax.set_extent(extent, geodetic)
>>> geoms = list(dataset.intersecting_geometries(extent))
>>>
>>> # Add the Stamen aerial imagery at zoom level 7.
>>>
>>> # Transform the geodetic coordinates of the coastlines into the two
>>> # projections of differing ellipses.
>>> wgs84_geoms = [geom_transform(transform_fn_factory(wgs84, geodetic),
>>>                               geom) for geom in geoms]
>>> sphere_geoms = [geom_transform(transform_fn_factory(sphere, geodetic),
>>>                                geom) for geom in geoms]
>>>
>>> # Using these differently referenced geometries, assume that they are
>>> # both referenced to WGS84.
>>>
>>> # Create a legend for the coastlines.
>>> legend_artists = [Line([0], [0], color=color, linewidth=3)
>>>                   for color in ('white', 'gray')]
>>> legend_texts = ['Correct ellipse\n(WGS84)', 'Incorrect ellipse\n(sphere)']
>>> legend = ax.legend(legend_artists, legend_texts, fancybox=True,
>>>                    loc='lower left', framealpha=0.75)
>>> legend.legendPatch.set_facecolor('wheat')
>>>
>>> # Create an inset GeoAxes showing the location of the Solomon Islands.
>>> sub_ax = fig.add_axes([0.7, 0.625, 0.2, 0.2],
>>>                       projection=ccrs.PlateCarree())
>>> sub_ax.set_extent([110, 180, -50, 10], geodetic)
>>>
>>> # Make a nice border around the inset axes.
>>> effect = Stroke(linewidth=4, foreground='wheat', alpha=0.5)
>>> sub_ax.outline_patch.set_path_effects([effect])
>>>
>>> # Add the land, coastlines and the extent of the Solomon Islands.
>>> sub_ax.coastlines()
>>> extent_box = sgeom.box(extent[0], extent[2], extent[1], extent[3])
>>>                       edgecolor='blue', linewidth=2)
>>>
>>> plt.show()

/usr/lib/python3/dist-packages/cartopy/io/__init__.py:260: DownloadWarning: Downloading: http://naciscdn.org/naturalearth/10m/physical/ne_10m_coastline.zip


## 11.5.2. Tissot’s indicatrix (蒂索氏指数)¶

>>> import matplotlib.pyplot as plt
>>>
>>> import cartopy.crs as ccrs
>>>
>>> fig = plt.figure(figsize=(10, 5))
>>> ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
>>>
>>> # make the map global rather than have it zoom in to
>>> # the extents of any plotted data
>>> ax.set_global()
>>>
>>> ax.stock_img()
>>> ax.coastlines()
>>>
>>> ax.tissot(facecolor='orange', alpha=0.4)
>>>
>>> plt.show()

/usr/lib/python3/dist-packages/cartopy/mpl/geoaxes.py:632: UserWarning: Approximating coordinate system <cartopy._crs.Geodetic object at 0x7f383397d518> with the PlateCarree projection.
'PlateCarree projection.'.format(crs))


## 11.5.3. UN 国旗¶

Matplotlib 与 Cartopy 的方位等距投影相结合，生成联合国国旗。

>>> import cartopy.crs as ccrs
>>> import cartopy.feature as cfeature
>>> import matplotlib.pyplot as plt
>>> from matplotlib.patches import PathPatch
>>> import matplotlib.path
>>> import matplotlib.ticker
>>> from matplotlib.transforms import BboxTransform, Bbox
>>> import numpy as np
>>>
>>>
>>> # When drawing the flag, we can either use white filled land, or be a little
>>> # more fancy and use the Natural Earth shaded relief imagery.
>>> filled_land = True
>>>
>>>
>>> def olive_path():
>>>     """
>>>     Return a Matplotlib path representing a single olive branch from the
>>>     UN Flag. The path coordinates were extracted from the SVG at
>>>     https://commons.wikimedia.org/wiki/File:Flag_of_the_United_Nations.svg.
>>>
>>>     """
>>>     olives_verts = np.array(
>>>         [[0,   2,   6,   9,  30,  55,  79,  94, 104, 117, 134, 157, 177,
>>>           188, 199, 207, 191, 167, 149, 129, 109,  87,  53,  22,   0, 663,
>>>           245, 223, 187, 158, 154, 150, 146, 149, 154, 158, 181, 184, 197,
>>>           181, 167, 153, 142, 129, 116, 119, 123, 127, 151, 178, 203, 220,
>>>           237, 245, 663, 280, 267, 232, 209, 205, 201, 196, 196, 201, 207,
>>>           211, 224, 219, 230, 220, 212, 207, 198, 195, 176, 197, 220, 239,
>>>           259, 277, 280, 663, 295, 293, 264, 250, 247, 244, 240, 240, 243,
>>>           244, 249, 251, 250, 248, 242, 245, 233, 236, 230, 228, 224, 222,
>>>           234, 249, 262, 275, 285, 291, 295, 296, 295, 663, 294, 293, 292,
>>>           289, 294, 277, 271, 269, 268, 265, 264, 264, 264, 272, 260, 248,
>>>           245, 243, 242, 240, 243, 245, 247, 252, 256, 259, 258, 257, 258,
>>>           267, 285, 290, 294, 297, 294, 663, 285, 285, 277, 266, 265, 265,
>>>           265, 277, 266, 268, 269, 269, 269, 268, 268, 267, 267, 264, 248,
>>>           235, 232, 229, 228, 229, 232, 236, 246, 266, 269, 271, 285, 285,
>>>           663, 252, 245, 238, 230, 246, 245, 250, 252, 255, 256, 256, 253,
>>>           249, 242, 231, 214, 208, 208, 227, 244, 252, 258, 262, 262, 261,
>>>           262, 264, 265, 252, 663, 185, 197, 206, 215, 223, 233, 242, 237,
>>>           237, 230, 220, 202, 185, 663],
>>>          [8,   5,   3,   0,  22,  46,  46,  46,  35,  27,  16,  10,  18,
>>>           22,  28,  38,  27,  26,  33,  41,  52,  52,  52,  30,   8, 595,
>>>           77,  52,  61,  54,  53,  52,  53,  55,  55,  57,  65,  90, 106,
>>>           96,  81,  68,  58,  54,  51,  50,  51,  50,  44,  34,  43,  48,
>>>           61,  77, 595, 135, 104, 102,  83,  79,  76,  74,  74,  79,  84,
>>>           90, 109, 135, 156, 145, 133, 121, 100,  77,  62,  69,  67,  80,
>>>           92, 113, 135, 595, 198, 171, 156, 134, 129, 124, 120, 123, 126,
>>>           129, 138, 149, 161, 175, 188, 202, 177, 144, 116, 110, 105,  99,
>>>           108, 116, 126, 136, 147, 162, 173, 186, 198, 595, 249, 255, 261,
>>>           267, 241, 222, 200, 192, 183, 175, 175, 175, 175, 199, 221, 240,
>>>           245, 250, 256, 245, 233, 222, 207, 194, 180, 172, 162, 153, 154,
>>>           171, 184, 202, 216, 233, 249, 595, 276, 296, 312, 327, 327, 327,
>>>           327, 308, 284, 262, 240, 240, 239, 239, 242, 244, 247, 265, 277,
>>>           290, 293, 296, 300, 291, 282, 274, 253, 236, 213, 235, 252, 276,
>>>           595, 342, 349, 355, 357, 346, 326, 309, 303, 297, 291, 290, 297,
>>>           304, 310, 321, 327, 343, 321, 305, 292, 286, 278, 270, 276, 281,
>>>           287, 306, 328, 342, 595, 379, 369, 355, 343, 333, 326, 318, 328,
>>>           340, 349, 366, 373, 379, 595]]).T
>>>     olives_codes = np.array([1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
>>>                              4, 4, 4, 4, 4, 4, 4, 4, 4, 79, 1, 4, 4, 4, 4, 4,
>>>                              4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
>>>                              4, 4, 4, 4, 4, 4, 79, 1, 4, 4, 4, 4, 4, 4, 2, 4,
>>>                              4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
>>>                              4, 79, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
>>>                              4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
>>>                              4, 79, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
>>>                              4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 2, 4,
>>>                              4, 4, 4, 4, 4, 79, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4,
>>>                              2, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
>>>                              4, 4, 4, 4, 4, 4, 79, 1, 4, 4, 4, 4, 4, 4, 4, 4,
>>>                              4, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
>>>                              4, 4, 4, 4, 79, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
>>>                              4, 4, 79], dtype=np.uint8)
>>>
>>>     return matplotlib.path.Path(olives_verts, olives_codes)
>>>
>>>
>>>
>>> blue = '#4b92db'
>>>
>>> # We're drawing a flag with a 3:5 aspect ratio.
>>> fig = plt.figure(figsize=[7.5, 4.5], facecolor=blue)
>>> # Put a blue background on the figure.
>>> blue_background = PathPatch(matplotlib.path.Path.unit_rectangle(),
>>>                             transform=fig.transFigure, color=blue,
>>>                             zorder=-1)
>>> fig.patches.append(blue_background)
>>>
>>> # Set up the Azimuthal Equidistant and Plate Carree projections
>>> # for later use.
>>> az_eq = ccrs.AzimuthalEquidistant(central_latitude=90)
>>> pc = ccrs.PlateCarree()
>>>
>>> # Pick a suitable location for the map (which is in an Azimuthal
>>> # Equidistant projection).
>>> ax = fig.add_axes([0.25, 0.24, 0.5, 0.54], projection=az_eq)
>>>
>>> # The background patch and outline patch are not needed in this example.
>>> ax.background_patch.set_facecolor('none')
>>> ax.outline_patch.set_edgecolor('none')
>>>
>>> # We want the map to go down to -60 degrees latitude.
>>> ax.set_extent([-180, 180, -60, 90], ccrs.PlateCarree())
>>>
>>> # Importantly, we want the axes to be circular at the -60 latitude
>>> # rather than cartopy's default behaviour of zooming in and becoming
>>> # square.
>>> _, patch_radius = az_eq.transform_point(0, -60, pc)
>>> ax.set_boundary(circular_path)
>>>
>>> if filled_land:
>>>         cfeature.LAND, facecolor='white', edgecolor='none')
>>> else:
>>>     ax.stock_img()
>>>
>>> gl = ax.gridlines(crs=pc, linewidth=3, color='white', linestyle='-')
>>> # Meridians every 45 degrees, and 5 parallels.
>>> gl.xlocator = matplotlib.ticker.FixedLocator(np.arange(0, 361, 45))
>>> parallels = np.linspace(-60, 70, 5, endpoint=True)
>>> gl.ylocator = matplotlib.ticker.FixedLocator(parallels)
>>>
>>> # Now add the olive branches around the axes. We do this in normalised
>>> # figure coordinates
>>> olive_leaf = olive_path()
>>>
>>> olives_bbox = Bbox.null()
>>> olives_bbox.update_from_path(olive_leaf)
>>>
>>> # The first olive branch goes from left to right.
>>> olive1_axes_bbox = Bbox([[0.45, 0.15], [0.725, 0.75]])
>>> olive1_trans = BboxTransform(olives_bbox, olive1_axes_bbox)
>>>
>>> # THe second olive branch goes from right to left (mirroring the first).
>>> olive2_axes_bbox = Bbox([[0.55, 0.15], [0.275, 0.75]])
>>> olive2_trans = BboxTransform(olives_bbox, olive2_axes_bbox)
>>>
>>> olive1 = PathPatch(olive_leaf, facecolor='white', edgecolor='none',
>>>                    transform=olive1_trans + fig.transFigure)
>>> olive2 = PathPatch(olive_leaf, facecolor='white', edgecolor='none',
>>>                    transform=olive2_trans + fig.transFigure)
>>>
>>> fig.patches.append(olive1)
>>> fig.patches.append(olive2)
>>>
>>> plt.show()
>>>
>>>


## 11.5.4. 在偏心椭圆上显示数据¶

>>> from urllib.request import urlopen
>>> from io import BytesIO
>>>
>>> import cartopy.crs as ccrs
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from PIL import Image
>>>
>>>
>>> def vesta_image():
>>>     """
>>>     Return an image of Vesta's topography.
>>>
>>>     Image credit: NASA/JPL-Caltech/UCLA/MPS/DLR/IDA/PSI
>>>
>>>     Returns
>>>     -------
>>>     img : numpy array
>>>         The pixels of the image in a numpy array.
>>>     img_proj : cartopy CRS
>>>         The rectangular coordinate system of the image.
>>>     img_extent : tuple of floats
>>>         The extent of the image (x0, y0, x1, y1) referenced in
>>>         the img_proj coordinate system.
>>>
>>>     """
>>>     url = 'https://www.nasa.gov/sites/default/files/pia17037.jpg'
>>>     raw_image = Image.open(img_handle)
>>>     # The image is extremely high-resolution, which takes a long time to
>>>     # plot. Sub-sampling reduces the time taken to plot while not
>>>     # significantly altering the integrity of the result.
>>>     smaller_image = raw_image.resize([raw_image.size[0] // 10,
>>>                                       raw_image.size[1] // 10])
>>>     img = np.asarray(smaller_image)
>>>     # We define the semimajor and semiminor axes, but must also tell the
>>>     # globe not to use the WGS84 ellipse, which is its default behaviour.
>>>     img_globe = ccrs.Globe(semimajor_axis=285000., semiminor_axis=229000.,
>>>                            ellipse=None)
>>>     img_proj = ccrs.PlateCarree(globe=img_globe)
>>>     img_extent = (-895353.906273091, 895353.906273091,
>>>                   447676.9531365455, -447676.9531365455)
>>>     return img, img_globe, img_proj, img_extent
>>>
>>>
>>>
>>> img, globe, crs, extent = vesta_image()
>>> projection = ccrs.Geostationary(globe=globe)
>>>
>>> fig = plt.figure()
>>> ax = fig.add_subplot(1, 1, 1, projection=projection)
>>> ax.imshow(img, transform=crs, extent=extent)
>>> fig.text(.075, .012, "Image credit: NASA/JPL-Caltech/UCLA/MPS/DLR/IDA/PSI",
>>>          bbox={'facecolor': 'w', 'edgecolor': 'k'})
>>> plt.show()


## 11.5.5. 显示UTM投影的所有60个区域¶

>>> import cartopy.crs as ccrs
>>> import matplotlib.pyplot as plt
>>>
>>>
>>>
>>> # Create a list of integers from 1 - 60
>>> zones = range(1, 61)
>>>
>>> # Create a figure
>>> fig = plt.figure(figsize=(18, 6))
>>>
>>> # Loop through each zone in the list
>>> for zone in zones:
>>>
>>>     # Add GeoAxes object with specific UTM zone projection to the figure
>>>         1, len(zones), zone,
>>>         projection=ccrs.UTM(
>>>             zone=zone,
>>>             southern_hemisphere=True
>>>         )
>>>     )
>>>
>>>     # Add coastlines, gridlines and zone number for the subplot
>>>     ax.coastlines(resolution='110m')
>>>     ax.gridlines()
>>>     ax.set_title(zone)
>>>
>>> # Add a supertitle for the figure
>>> fig.suptitle("UTM Projection - Zones")
>>>
>>> # Display the figure
>>> plt.show()