标签：30 ticks add 分布图 小白 ccrs cartopy ax ds
小白学习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，没事继续跑着
测试所需数据可以在评论区留下网盘地址，作者给你们发，希望看过如果觉得有用的朋友点个赞评论一下，感谢哟。
标签：30,ticks,add,分布图,小白,ccrs,cartopy,ax,ds 来源： https://blog.csdn.net/weixin_42372313/article/details/113665724
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。