小白学习cartopy气象画地图的第二天（中国区域，陆地温度分布图） |
您所在的位置：网站首页 › 如何画平面区域图 › 小白学习cartopy气象画地图的第二天（中国区域，陆地温度分布图） |
小白学习cartopy气象画地图的第二天（中国区域，陆地温度分布图）
|
小白学习cartopy气象画地图的第二天（中国区域，陆地温度分布图）
首先，还是先感谢一下公众号（气象学家），因为代码和测试数据都是他的，我只是拿过来，然后慢慢读懂它，最后再以小白的角度把每一行代码解释给大家听，拾人牙慧而已。话不多说，开始！ 这次用到的就比较多啦 第二步先写一些功能函数 #给画刻度用来辅助，确定四周 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.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) & (lon10) & (lat'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.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) & (lon10) & (lat |
【本文地址】
公司简介
联系我们
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |