小白学习cartopy气象画地图的第二天（中国区域，陆地温度分布图）
首先，还是先感谢一下公众号（气象学家），因为代码和测试数据都是他的，我只是拿过来，然后慢慢读懂它，最后再以小白的角度把每一行代码解释给大家听，拾人牙慧而已。话不多说，开始！
第一步导入库
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from copy import copy
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
import shapely.geometry as sgeom
这次用到的就比较多啦
第二步先写一些功能函数
#给画刻度用来辅助，确定四周
def find_side(ls, side):
"""
Given a shapely LineString which is assumed to be rectangular, return the
line corresponding to a given side of the rectangle.
"""
minx, miny, maxx, maxy = ls.bounds
points = {'left': [(minx, miny), (minx, maxy)],
'right': [(maxx, miny), (maxx, maxy)],
'bottom': [(minx, miny), (maxx, miny)],
'top': [(minx, maxy), (maxx, maxy)],}
return sgeom.LineString(points[side])
#用来画坐标轴的刻度（包括经度和纬度）
def _lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor):
"""得到一个兰伯特正形投影的轴的刻度位置和标签."""
outline_patch = sgeom.LineString(ax.outline_patch.get_path().vertices.tolist())
axis = find_side(outline_patch, tick_location)
n_steps = 30
extent = ax.get_extent(ccrs.PlateCarree())
_ticks = []
for t in ticks:
xy = line_constructor(t, n_steps, extent)
proj_xyz = ax.projection.transform_points(ccrs.Geodetic(), xy[:, 0], xy[:, 1])
xyt = proj_xyz[..., :2]
ls = sgeom.LineString(xyt.tolist())
locs = axis.intersection(ls)
if not locs:
tick = [None]
else:
tick = tick_extractor(locs.xy)
_ticks.append(tick[0])
# Remove ticks that aren't visible:
ticklabels = copy(ticks)
while True:
try:
index = _ticks.index(None)
except ValueError:
break
_ticks.pop(index)
ticklabels.pop(index)
return _ticks, ticklabels
#设置x轴标签（用来画纬度）
def lambert_xticks(ax, ticks):
"""Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
te = lambda xy: xy[0]
lc = lambda t, n, b: np.vstack((np.zeros(n) + t, np.linspace(b[2], b[3], n))).T
xticks, xticklabels = _lambert_ticks(ax, ticks, 'bottom', lc, te)
ax.xaxis.tick_bottom()
ax.set_xticks(xticks)
ax.set_xticklabels([ax.xaxis.get_major_formatter()(xtick) for xtick in xticklabels])
#设置y轴标签（用来画经度）
def lambert_yticks(ax, ticks):
"""Draw ricks on the left y-axis of a Lamber Conformal projection."""
te = lambda xy: xy[1]
lc = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n) + t)).T
yticks, yticklabels = _lambert_ticks(ax, ticks, 'left', lc, te)
ax.yaxis.tick_left()
ax.set_yticks(yticks)
ax.set_yticklabels([ax.yaxis.get_major_formatter()(ytick) for ytick in yticklabels])
#掩膜函数
def mask(ds, label='land'):
landsea = xr.open_dataset('landsea.nc')
landsea = landsea['LSMASK']
#ds和地形数据分辨率不一致，需将地形数据插值
landsea = landsea.interp(lat=ds.latitude.values, lon=ds.longitude.values)
#利用地形掩盖海陆数据
ds.coords['mask'] = (('latitude', 'longitude'), landsea.values)
# print(ds.mask)
if label == 'land':
ds = ds.where(ds.mask < 0.8) #可以尝试调整default：0.8
elif label == 'ocean':
ds = ds.where(ds.mask > 0.3) #可以尝试调整default：0.2
return ds
前四个是用来画坐标轴刻度的，也就是经纬度刻度，纬度是Y轴，经度是X轴。第五个mask()是用来做掩膜的，其实就是剔除不想要的数据，这里是用来剔除海洋数据的。
第三步画地图底图
这一步在我上一文章里面已经讲过了，不再赘述了，但是里面有一些细节上的不同我已经通过注释的方式写出来了，慢慢品
#读取CN-border-La.dat文件
with open('CN-border-La.dat') as src:
context = src.read()
blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]
#设置画图各种参数
fig = plt.figure(figsize=[10, 8],frameon=True)
ax = fig.add_axes([0.08, 0.05, 0.8, 0.94], projection=ccrs.LambertConformal(central_latitude=90, central_longitude=105))
# ccrs.LambertConformal( central_longitude=-96.0, central_latitude=39.0,
# false_easting=0.0, false_northing=0.0,
# secant_latitudes=None, standard_parallels=None,
# globe=None, cutoff=-30)
#参数含义：
# central_longitude:*经度。默认为-96。
# central_latitude:*纬度。默认为39。
# false_easting:平面原点的X偏移量，单位为米。默认值为0。
# false_northing:Y距平面原点的距离，单位为米。默认值为0。
# secant_latitudes:正割纬度。该关键字在v0.12中已直接弃用， 被‘standard parallels（标准平行）’取代。默认为没有。
# standard_parallels:标准并行纬度(s)。默认为(33,45)。
# globe:类:“cartopy.crs.Globe”。如果省略，则为默认创建全局变量。
# cutoff:可选，地图截止纬度。地图一直延伸到*极点对面的无穷远处，所以我们必须在那之前把地图剪掉。该值为0将绘制半个地球。默认为-30。
# 画海，陆地，河流，湖泊
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
# 画国界
for line in borders:
ax.plot(line[0::2], line[1::2], '-', lw=1.5, color='k',
transform=ccrs.Geodetic())
# 框出区域
ax.set_extent([80, 130, 15, 55],crs=ccrs.PlateCarree())
#必须调用.draw()函数，用来画刻度的边框，没有这条，下面会报错，没法画刻度
fig.canvas.draw()
#设置刻度值，画经纬网格
xticks = [55, 65, 75, 85, 95, 105, 115, 125, 135, 145, 155, 165]
yticks = [0 , 5 , 10, 15, 20, 25 , 30 , 35 , 40 , 45 , 50 , 55 , 60 , 65]
ax.gridlines(xlocs=xticks, ylocs=yticks,linestyle='--',lw=1,color='dimgrey')
# 设置经纬网格的端点（也是用来配合画刻度的，注释掉以后刻度不能正常显示了）
ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)
#画经纬度刻度
lambert_xticks(ax, xticks)
lambert_yticks(ax, yticks)
#画南海
sub_ax = fig.add_axes([0.592, 0.189, 0.14, 0.155],
projection=ccrs.LambertConformal(central_latitude=90,
central_longitude=115))
sub_ax.add_feature(cfeature.OCEAN.with_scale('50m'))
sub_ax.add_feature(cfeature.LAND.with_scale('50m'))
sub_ax.add_feature(cfeature.RIVERS.with_scale('50m'))
sub_ax.add_feature(cfeature.LAKES.with_scale('50m'))
for line in borders:
sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',
transform=ccrs.Geodetic())
sub_ax.set_extent([105, 125, 0, 25],crs=ccrs.PlateCarree())
友情提供：
ccrs.LambertConformal(central_longitude=-96.0,central_latitude=39.0,false_easting=0.0, false_northing=0.0, secant_latitudes=None, standard_parallels=None, globe=None, cutoff=-30)
参数含义：
central_longitude:*经度。默认为-96。
central_latitude:*纬度。默认为39。
false_easting:平面原点的X偏移量，单位为米。默认值为0。
false_northing:Y距平面原点的距离，单位为米。默认值为0。
secant_latitudes:正割纬度。该关键字在v0.12中已直接弃用， 被‘standard parallels（标准平行）’取代。默认为没有。
standard_parallels:标准并行纬度(s)。默认为(33,45)。
globe:类:“cartopy.crs.Globe”。如果省略，则为默认创建全局变量。
cutoff:地图截止纬度。地图一直延伸到*极点对面的无穷远处，所以我们必须在那之前把地图剪掉。该值为0将绘制半个地球。默认为-30。
第四步读取温度数据（NC格式）
#加载NC资料，用相对路径可能会报错，如果错了可以试试绝对路径
#ds指Dataset，另外还有da指DataArray，是xarray的一种数据结构
ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
#拿到数据基本信息，经纬度，时间，温度
lat = ds.latitude
lon = ds.longitude
time = ds.time
temp = (ds['t2m'] - 273.15) # 把温度转换为℃
#设置经纬度范围
lon_range = lon[(lon>60) & (lon<150)]
lat_range = lat[(lat>10) & (lat<60)]
#.sel()的意思就是通过取出2018-02-01那天，我们需要的经纬度区域的数据
temp_region = temp.sel(longitude=lon_range, latitude=lat_range, time='2018-02-01')
详细的每一行是啥意思我也已经写在注释里面了
第五步画图
设置掩膜
#把海洋部分用掩膜覆盖，其实就是用mask函数剔除海洋部分的数据
temp_mask = mask(temp_region, 'ocean')
设置色标
cbar_kwargs = {
'orientation': 'vertical',
'label': '2m temperature (℃)',
'shrink': 0.8,
'ticks': np.arange(-30, 30 + 5, 5),
'pad': 0.05,
'shrink': 0.65
}
设置数据精度
#设置画数据的精度，这里设置了从-30到30度，按1度画温度
levels = np.arange(-30, 30 + 1, 1)
画图
在这里插入代码片
设置色标
#本质上是一个画等值线并填色的函数，被cartopy封装了一下，参数的含义可以去查matplotlib.pyplot.contourf。
temp_mask.plot.contourf(ax=ax, levels=levels, cmap='Spectral_r',
cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())
#参数含义：levels设置等值线密度，cmap颜色，cbar_kwargs色标，transform变换
设置标题
# 设置标题的在代码中放置的位置很关键，注意不要放置在小图上或者新建画框了。
ax.set_title('China Map 2m Temperature',color='blue',fontsize= 20)
最后完整代码
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from copy import copy
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
import shapely.geometry as sgeom
#给画刻度用来辅助，确定四周
def find_side(ls, side):
"""
Given a shapely LineString which is assumed to be rectangular, return the
line corresponding to a given side of the rectangle.
"""
minx, miny, maxx, maxy = ls.bounds
points = {'left': [(minx, miny), (minx, maxy)],
'right': [(maxx, miny), (maxx, maxy)],
'bottom': [(minx, miny), (maxx, miny)],
'top': [(minx, maxy), (maxx, maxy)],}
return sgeom.LineString(points[side])
#用来画坐标轴的刻度（包括经度和纬度）
def _lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor):
"""得到一个兰伯特正形投影的轴的刻度位置和标签."""
outline_patch = sgeom.LineString(ax.outline_patch.get_path().vertices.tolist())
axis = find_side(outline_patch, tick_location)
n_steps = 30
extent = ax.get_extent(ccrs.PlateCarree())
_ticks = []
for t in ticks:
xy = line_constructor(t, n_steps, extent)
proj_xyz = ax.projection.transform_points(ccrs.Geodetic(), xy[:, 0], xy[:, 1])
xyt = proj_xyz[..., :2]
ls = sgeom.LineString(xyt.tolist())
locs = axis.intersection(ls)
if not locs:
tick = [None]
else:
tick = tick_extractor(locs.xy)
_ticks.append(tick[0])
# Remove ticks that aren't visible:
ticklabels = copy(ticks)
while True:
try:
index = _ticks.index(None)
except ValueError:
break
_ticks.pop(index)
ticklabels.pop(index)
return _ticks, ticklabels
#设置x轴标签（用来画纬度）
def lambert_xticks(ax, ticks):
"""Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
te = lambda xy: xy[0]
lc = lambda t, n, b: np.vstack((np.zeros(n) + t, np.linspace(b[2], b[3], n))).T
xticks, xticklabels = _lambert_ticks(ax, ticks, 'bottom', lc, te)
ax.xaxis.tick_bottom()
ax.set_xticks(xticks)
ax.set_xticklabels([ax.xaxis.get_major_formatter()(xtick) for xtick in xticklabels])
#设置y轴标签（用来画经度）
def lambert_yticks(ax, ticks):
"""Draw ricks on the left y-axis of a Lamber Conformal projection."""
te = lambda xy: xy[1]
lc = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n) + t)).T
yticks, yticklabels = _lambert_ticks(ax, ticks, 'left', lc, te)
ax.yaxis.tick_left()
ax.set_yticks(yticks)
ax.set_yticklabels([ax.yaxis.get_major_formatter()(ytick) for ytick in yticklabels])
#掩膜函数
def mask(ds, label='land'):
landsea = xr.open_dataset('landsea.nc')
landsea = landsea['LSMASK']
#ds和地形数据分辨率不一致，需将地形数据插值
landsea = landsea.interp(lat=ds.latitude.values, lon=ds.longitude.values)
#利用地形掩盖海陆数据
ds.coords['mask'] = (('latitude', 'longitude'), landsea.values)
# print(ds.mask)
if label == 'land':
ds = ds.where(ds.mask < 0.8) #可以尝试调整default：0.8
elif label == 'ocean':
ds = ds.where(ds.mask > 0.3) #可以尝试调整default：0.2
return ds
#----------------------分割线（上面是功能函数，后面会用到）------------------------
#读取CN-border-La.dat文件
with open('CN-border-La.dat') as src:
context = src.read()
blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]
#设置画图各种参数
fig = plt.figure(figsize=[10, 8],frameon=True)
ax = fig.add_axes([0.08, 0.05, 0.8, 0.94], projection=ccrs.LambertConformal(central_latitude=90, central_longitude=105))
# ccrs.LambertConformal( central_longitude=-96.0, central_latitude=39.0,
# false_easting=0.0, false_northing=0.0,
# secant_latitudes=None, standard_parallels=None,
# globe=None, cutoff=-30)
# 参数含义：
# central_longitude:*经度。默认为-96。
# central_latitude:*纬度。默认为39。
# false_easting:平面原点的X偏移量，单位为米。默认值为0。
# false_northing:Y距平面原点的距离，单位为米。默认值为0。
# secant_latitudes:正割纬度。该关键字在v0.12中已直接弃用， 被‘standard parallels（标准平行）’取代。默认为没有。
# standard_parallels:标准并行纬度(s)。默认为(33,45)。
# globe:类:“cartopy.crs.Globe”。如果省略，则为默认创建全局变量。
# cutoff:可选，地图截止纬度。地图一直延伸到*极点对面的无穷远处，所以我们必须在那之前把地图剪掉。该值为0将绘制半个地球。默认为-30。
# 画海，陆地，河流，湖泊
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
# 画国界
for line in borders:
ax.plot(line[0::2], line[1::2], '-', lw=1.5, color='k',
transform=ccrs.Geodetic())
# 框出区域
ax.set_extent([80, 130, 15, 55],crs=ccrs.PlateCarree())
#必须调用.draw()函数，用来画刻度的边框，没有这条，下面会报错，没法画刻度
fig.canvas.draw()
#设置刻度值，画经纬网格
xticks = [55, 65, 75, 85, 95, 105, 115, 125, 135, 145, 155, 165]
yticks = [0 , 5 , 10, 15, 20, 25 , 30 , 35 , 40 , 45 , 50 , 55 , 60 , 65]
ax.gridlines(xlocs=xticks, ylocs=yticks,linestyle='--',lw=1,color='dimgrey')
# 设置经纬网格的端点（也是用来配合画刻度的，注释掉以后刻度不能正常显示了）
ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)
#画经纬度刻度
lambert_xticks(ax, xticks)
lambert_yticks(ax, yticks)
#画南海
sub_ax = fig.add_axes([0.592, 0.189, 0.14, 0.155],
projection=ccrs.LambertConformal(central_latitude=90,
central_longitude=115))
sub_ax.add_feature(cfeature.OCEAN.with_scale('50m'))
sub_ax.add_feature(cfeature.LAND.with_scale('50m'))
sub_ax.add_feature(cfeature.RIVERS.with_scale('50m'))
sub_ax.add_feature(cfeature.LAKES.with_scale('50m'))
for line in borders:
sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',
transform=ccrs.Geodetic())
sub_ax.set_extent([105, 125, 0, 25],crs=ccrs.PlateCarree())
#----------------------分割线（至此，地图地图加载完毕）------------------------
#加载NC资料，用相对路径可能会报错，如果错了可以试试绝对路径
#ds指Dataset，另外还有da指DataArray，是xarray的一种数据结构
ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
#拿到数据基本信息，经纬度，时间，温度
lat = ds.latitude
lon = ds.longitude
time = ds.time
temp = (ds['t2m'] - 273.15) # 把温度转换为℃
#设置经纬度范围
lon_range = lon[(lon>60) & (lon<150)]
lat_range = lat[(lat>10) & (lat<60)]
#.sel()的意思就是通过取出2018-02-01那天，我们需要的经纬度区域的数据
temp_region = temp.sel(longitude=lon_range, latitude=lat_range, time='2018-02-01')
#把海洋部分用掩膜覆盖，其实就是用mask函数剔除海洋部分的数据
temp_mask = mask(temp_region, 'ocean')
#设置色标的
cbar_kwargs = {
'orientation': 'vertical',
'label': '2m temperature (℃)',
'shrink': 0.8,
'ticks': np.arange(-30, 30 + 5, 5),
'pad': 0.05,
'shrink': 0.65
}
#设置画数据的精度，这里设置了从-30到30度，按1度画温度
levels = np.arange(-30, 30 + 1, 1)
#本质上是一个画等值线并填色的函数，被cartopy封装了一下，参数的含义可以去查matplotlib.pyplot.contourf。
temp_mask.plot.contourf(ax=ax, levels=levels, cmap='Spectral_r',
cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())
#参数含义：levels设置等值线密度，cmap颜色，cbar_kwargs色标，transform变换
#最后画个标题
# 设置标题的在代码中放置的位置很关键，注意不要放置在小图上或者新建画框了。
ax.set_title('China Map 2m Temperature',color='blue',fontsize= 20)
#----------------------分割线（至此数据也画完了）------------------------
#出图
plt.show()
代码可能会跑一个警告DeprecationWarning，没事继续跑着
测试所需数据可以在评论区留下网盘地址，作者给你们发，希望看过如果觉得有用的朋友点个赞评论一下，感谢哟。