摘要: 今天我们将学习 Cartopy,它是在 python 中绘制地图的最常见的包之一。一个受欢迎且功能强大的库是 Basemap。Basemap 即将消失,并在不久的将来被 Cartopy 取代,因此,建议您花时间来学习。 Cartopy 利用强大的PROJ.4,...
今天我们将学习 Cartopy,它是在 python 中绘制地图的最常见的包之一。一个受欢迎且功能强大的库是 Basemap。Basemap 即将消失,并在不久的将来被 Cartopy 取代,因此,建议您花时间来学习。
Cartopy 利用强大的PROJ.4,numpy和shapely库,并包括基于Matplotlib构建的编程接口,用于创建出版质量地图。cartopy 的主要特点是其面向对象的投影定义,以及在这些投影之间转换点、线、向量、多边形和图像的能力。
Cartopy Projections 和其他参考系统
在 Cartopy 中,每个投影都是一个类。尽管 Cartopy 对合理的默认设置持固执己见的立场,但大多数类型的投影都能以特定于投影的方式进行配置,以创建一个 Plate Carree 投影实例,需要 cartopy 的 crs 模块。这通常作为 ccrs(Cartopy 坐标参考系统)导入。
在导入绘制地图所需的模块之前,首先要使用reticulate包在 Rstudio 中初始化 python 的链接,并设置 seesion 获取 python 函数和包的环境,须在 R 块中加载网状包并设置环境;
require(reticulate)
## Loading required package: reticulate
use_python("c:/Python/Anaconda3/")
利用import导入一些Python模块。确保插入Pyhon块以加载这些模块。
import numpy as np
import pandas as pd
import cartopy.crs as ccrs
import cartopy
import matplotlib.pyplot as plt
使用下面的代码块访问Mollweide投影;
ccrs.Mollweide()
## <cartopy.crs.Mollweide object at 0x0000000028A475E8>
绘制地图
Cartopy可选地依赖于matplotlib,每个投影都知道如何创建能够表示其自身的matplotlib Axes(或AxesSubploy)。
投影创建的轴是cartopy.mpl.Geoaxes.GeoAx。 Axes 子类覆盖了matplotlib的一些现有方法,并添加了许多非常有用的方法来绘制地图。
稍后将很快返回并查看这些方法,但首先,让我们看到cartopy+matplotlib 的实际效果:
import matplotlib.pyplot as plt plt.axes(projection=ccrs.PlateCarree()) plt.show()

令人印象深刻的是,创建的轴确实是GeoAx[SubPlot]实例之一。在标准matplotlib axes类之上添加的最有用的方法之一是CoastLine方法。在没有任何争议的情况下,它将向地图中添加自然地球1:110000000比例尺的海岸线数据。
plt.figure() ax = plt.axes(projection=ccrs.PlateCarree()) ax.coastlines() plt.show()

同样可以使用现有的许多方法之一创建 matplotlib 子图。 例如,可以使用 plt.subplots 函数:
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()}) ax.coastlines() plt.show
## <function make_python_function.<locals>.python_function at 0x0000000026489438>
投影类有一些选项可以用来自定义地图,以便非洲位于中心位置
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=0)) ax.coastlines() plt.show()

GeoAx的实用方法
Cartopy.mpl.Geoaxes.GeoAx类添加了许多有用的方法。让我们来看一下:
Ⅰ. set_global
-尽可能缩小地图
Ⅱ. set_extent
-将地图缩放到给定的边界框
Ⅲ. gridlines
-向轴添加经纬线(也可以添加标注)
Ⅳ. coastlines
-将自然地球海岸线添加到轴线
Ⅴ. stock_img
-将低分辨率的自然地球背景图像添加到轴
Ⅵ. imshow
-将图像(NumPy数组)添加到轴
Ⅶ. add_geometries
-将几何集合(形状)添加到轴
更多不同全球预测的示例
projections = [ccrs.PlateCarree(), ccrs.Robinson(), ccrs.Mercator(), ccrs.Orthographic(), ccrs.InterruptedGoodeHomolosine() ] for proj in projections: plt.figure() ax = plt.axes(projection=proj) ax.stock_img() ax.coastlines() ax.set_title(f'{type(proj)}') plt.show()





区域地图
为了创建区域地图,使用GeoAxis的Set_Extent方法来限制区域的大小。
central_lon =43, central_lat = -8.5 extent = [35, 50, -16.5, 0] ax = plt.axes(projection=ccrs.PlateCarree(central_lon, central_lat)) ax.set_extent(extent) ax.gridlines() ax.coastlines(resolution='50m') plt.show()
向地图添加要素
为了赋予地图更多的样式和细节,添加了cartopy.Feature对象。内置了许多有用的功能。这些“默认功能”的分辨率粗略计为 (110米)。
Ⅰ. cartopy.feature.BORDERS
国家/地区边界
Ⅱ. cartopy.feature.COASTLINE
海岸线,包括主要岛屿的 cartopy.feature.LAKES
天然和人工湖
Ⅲ. cartopy.feature.LAND
陆地多边形,包括主要岛屿
Ⅳ. cartopy.feature.OCEAN
海洋多边形
Ⅴ. cartopy.feature.RIVERS
单线排水系统,包括湖泊中心线
Ⅵ. cartopy.feature.STATES
(这种规模仅限于美国)
import cartopy.feature as cfeature import numpy as np central_lat = 37.5 central_lon = -96 extent = [28, 45, -25, 2] central_lon = np.mean(extent[:2]) central_lat = np.mean(extent[2:]) plt.figure(figsize=(6, 6)) ax = plt.axes(projection=ccrs.EquidistantConic(central_lon, central_lat)) ax.set_extent(extent) ax.add_feature(cartopy.feature.OCEAN) ax.add_feature(cartopy.feature.LAND, edgecolor='black') ax.add_feature(cartopy.feature.LAKES, edgecolor='black') ax.add_feature(cartopy.feature.RIVERS) ax.gridlines()
## <cartopy.mpl.gridliner.Gridliner object at 0x000000003791BD88>
plt.show()

rivers_50m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '10m') central_lat = 37.5 central_lon = -96 extent = [28, 45, -25, 2] central_lon = np.mean(extent[:2]) central_lat = np.mean(extent[2:]) plt.figure(figsize=(6, 6)) ax = plt.axes(projection=ccrs.EquidistantConic(central_lon, central_lat)) ax.set_extent(extent) ax.add_feature(cartopy.feature.OCEAN) ax.add_feature(cartopy.feature.LAND, edgecolor='black') ax.add_feature(cartopy.feature.LAKES, edgecolor='black') ax.add_feature(rivers_50m, facecolor='None', edgecolor='b') ax.gridlines()
## <cartopy.mpl.gridliner.Gridliner object at 0x0000000037A3F148>
plt.show()

绘制二维(栅格)数据
同样的原则也适用于2D数据。下面创建一些以正则经度/经度坐标定义的示例数据。在本例中,我们将加载全球海洋表面温度数据。数据以netcdf格式存储,因此需要加载:
import netCDF4 as nc
从netCDF4模块中调用一个数据集函数来读取文件
sst = nc.Dataset("e:/MatlabWorking/GHRSST/20150101.nc")
需要提取存储在文件中的不同变量。但在提取之前,必须查看文件的内部结构,并使用正确的名称标识变量。可以使用nc.Variables函数来实现这一点
sst.variables
我们注意到该文件是lon. lat, analysed_sst, 和time. 时间是单个间隔。提取变量,如下所示;
Time=sst.Variables[‘time’] LON=sst.Variables[‘lon’] Lat=sst.Variables[‘lat’] DATA=sst.Variables[‘Analysted_sst’]
由于数据在矩形网格中,所以还需要使用 np.meshgrid()
将 lon 和 lat 转换为矩形网格。
meshgrid 的目的是从 x 值数组和 y 值数组创建一个矩形网格。
lon2d, lat2d = np.meshgrid(lon, lat)
因为温度是以开尔文标度记录的,所以可以通过简单地减去 273 来简单地转换为度数
datar = data[0,:,:]-273 # datar = np.flipud(datar)
绘制全球海洋表面温度的空间分布图,如图1所示
plt.figure(figsize=(6,5)) ax = plt.axes(projection=ccrs.PlateCarree()) ax.set_global() # ax.set_extent([-170,170,-30,30]) ax.coastlines() ax.contourf(lon2d, lat2d, datar, cmap = "jet")
## <matplotlib.contour.QuadContourSet object at 0x0000000012567588>
ax.gridlines(draw_labels=True, xlocs=np.arange(-180,180,30))
## <cartopy.mpl.gridliner.Gridliner object at 0x00000000125DDC48>
plt.show()

图1:海面温度
import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from matplotlib.offsetbox import AnchoredText fig = plt.figure() # ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) ax = plt.axes(projection=ccrs.PlateCarree()) ax.set_extent([80, 170, -45, 30], crs=ccrs.PlateCarree()) # Put a background image on for nice sea rendering. ax.stock_img() # Create a feature for States/Admin 1 regions at 1:50m from Natural Earth states_provinces = cfeature.NaturalEarthFeature( category='cultural', name='admin_1_states_provinces_lines', scale='50m', facecolor='none') SOURCE = 'Natural Earth' LICENSE = 'public domain' ax.add_feature(cfeature.LAND) ax.add_feature(cfeature.COASTLINE) ax.add_feature(states_provinces, edgecolor='gray') # Add a text annotation for the license information to the # the bottom right corner. text = AnchoredText(r'$\mathcircled{{c}}$ {}; license: {}' ''.format(SOURCE, LICENSE), loc=4, prop={'size': 12}, frameon=True) ax.add_artist(text) plt.show()

fig = plt.figure() ax = plt.axes(projection=ccrs.PlateCarree()) ax.set_extent([20, 70, -40, 20], crs=ccrs.PlateCarree()) ax.add_feature(cfeature.LAND, edgecolor='gray') ax.add_feature(cfeature.BORDERS) ax.add_feature(cfeature.COASTLINE) ax.gridlines(xlocs=np.arange(20,70,10), draw_labels=True, crs=ccrs.PlateCarree())
## <cartopy.mpl.gridliner.Gridliner object at 0x0000000037A54748>
plt.show()

fig = plt.figure() ax = plt.axes(projection=ccrs.PlateCarree()) ax.set_extent([20, 70, -40, 20], crs=ccrs.PlateCarree()) # Put a background image on for nice sea rendering. ax.stock_img() ax.add_feature(cfeature.BORDERS) ax.add_feature(cfeature.COASTLINE) ax.gridlines(xlocs=np.arange(20,70,10), draw_labels=True, crs=ccrs.PlateCarree())
## <cartopy.mpl.gridliner.Gridliner object at 0x000000003815EC08>
plt.show()

fig = plt.figure() ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) ax.set_extent([-20, 60, -40, 45], crs=ccrs.PlateCarree()) ax.stock_img() ax.add_feature(cfeature.LAND) ax.add_feature(cfeature.OCEAN) ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle='dotted') ax.add_feature(cfeature.LAKES, alpha=0.5) ax.add_feature(cfeature.RIVERS) plt.show()

全球地图
一个简单地图的示例,用于比较两个位置之间的 Geodetic 和 Plate Carree 线。
fig = plt.figure(figsize=(8, 4)) ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson(central_longitude=0)) ax.set_global() ax.stock_img() ax.coastlines() ax.gridlines() # ax.plot(-0.08, 51.53, 'o', transform=ccrs.PlateCarree()) # ax.plot([-0.08, 132], [51.53, 43.17], transform=ccrs.PlateCarree()) # ax.plot([-0.08, 132], [51.53, 43.17], transform=ccrs.Geodetic())
## <cartopy.mpl.gridliner.Gridliner object at 0x000000003BD9A748>
plt.show()

填充轮廓
使用网格海洋表面温度,因为想看到太平洋和印度洋,所以使用了Mollweide投影,并将经度居中到130。
## read the netcdf file with Dataset function from netCDF4 module sst = nc.Dataset("e:/MatlabWorking/GHRSST/20050101.nc") ## sst.variables ## extract individual varibales in the nc file time = sst.variables['time'][:] lon = sst.variables['lon'][:] lat = sst.variables['lat'][:] sst = sst.variables['analysed_sst'] ## extract an array of the first day and convert temperature from Kelvin to Degree Celsius sstr = sst[0,:,:]-273 ## create a rectangular of the long and lat lon2d, lat2d = np.meshgrid(lon, lat) plt.figure(figsize=(8,4)) ax = plt.axes(projection=ccrs.Mollweide(central_longitude=160),) ax.set_global() # ax.set_extent([-170,170,-30,30]) ax.coastlines() ax.contourf(lon2d, lat2d, datar, transform=ccrs.PlateCarree(), cmap='nipy_spectral') # ax.gridlines(draw_labels=True, xlocs=np.arange(-180,180,30))# unsupported with Mollweide
## <matplotlib.contour.QuadContourSet object at 0x000000003BDEB308>
plt.show()

读取和绘制平均海平面异常图,如下图所示
## read mean sea level anomaly data msla = nc.Dataset("e:/MatlabWorking/Altimetry/msla_h/indian_ocean-twosat-msla-h_010193_311295.nc") ## msla.variables time = msla.variables['time'][:] lon = msla.variables['lon'][:] lat = msla.variables['lat'][:] sla = msla.variables['sla'] sla30 = sla[30,:,:] lon2d, lat2d = np.meshgrid(lon, lat) plt.figure(figsize=(6,5)) ax = plt.axes(projection=ccrs.Mollweide(central_longitude=60)) # ax.set_global() ax.set_extent([20,120,-50,20]) ax.coastlines() # ax.add_feature(cfeature.BORDERS, linestyle = "dotted") # uncomment to plot country boundaries ax.contourf(lon2d, lat2d, sla30, transform=ccrs.PlateCarree(), cmap='nipy_spectral')
## <matplotlib.contour.QuadContourSet object at 0x0000000037A6D548>
ax.gridlines(draw_labels=False, xlocs=np.arange(-180,180,30))
## <cartopy.mpl.gridliner.Gridliner object at 0x00000000349F7908>
plt.show()
