气象编程 | 数据筛选2-格点数据筛选

格点数据筛选

上一章节我们主要对站点数据进行了操作，本章节继续介绍对格点数据的筛选操作。格点数据有时候也叫栅格数据。气象上主要是再分析资料使用这种格式，主要有nc、grib、tif等格式。我们主要使用rioxarray来处理这些数据。

前期准备

导入一些常见的站点数据处理与可视化方法库包，基本设置等，这里我们选用了一份欧洲中心ERA5再分析天气资料：

import pandas as pd
import numpy as np
import xarray as xr
import rioxarray as rioxr
import geopandas as gpd
import cartopy.crs as ccrs
import cartopy.io.shapereader as cis
import cartopy.mpl.ticker as cmt
import cartopy.mpl.patch as cmp
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.colors as mcolors
import matplotlib.path as mpath
plt.rcParams['font.sans-serif']=['FangSong']
plt.rcParams['axes.unicode_minus']=False
data_path=r'C:\Users\Administrator\Desktop\20220727.grib'
ds=xr.open_dataset(data_path)

为何使用rioxarray

为什么我们不使用salem、regionmask等方法呢。因为这些方法或者安装麻烦，或者方法脱离xarray太远，而rioxarray是完全由xarray支持的，优势实在太大。 我们使用xarray读取资料后，只须在后面加上.rio就可以实现切换，ds还是xarray的Dataset，调用rio后就变成rioxarray的RasterDataset：

给再分析资料写入投影

为了进行后续的裁剪工作，我们需要给再分析资料写入投影，写入投影后，ds将会多出一个描述投影的spatial_ref。
为什么多写入epsg4326这个投影呢，因为绝大多数情况数据都是以这种投影存储的，他的中文名叫矩形投影,圆柱投影，在cartopy中叫PlateCarree。这种投影的经度为-180～180°，纬度为-90～90°。
在实际工作中遇到其他特殊情况进行特殊处理。

ds.rio.write_crs('EPSG:4326')

获取我们需要的地区的geometry

当shp文件仅含有一个区划时

这里我以湖北省为例进行筛选，首先读取湖北省shp文件：

hubei=gpd.read_file(path+r'\hubei\hubei.shp')

读取后可发现文件里只有湖北的几何数据。我们需要知道这份数据的投影：

hubei.crs

然后就可以裁剪数据了，rio.clip函数有两个必须的位置参数，第一个是用来裁剪的几何图形，在这里是湖北省的数据，第二个就是投影，这里使用该文件自带的投影,裁剪完成后，仍然输出xarray.Dataset格式，对后续绘图计算都非常友好，不在湖北省的数据被填充为nan，同时时间维度、气压层各高度都同时被裁剪：

ds_clip=ds.rio.write_crs('EPSG:4326').rio.clip(hubei.geometry.values,hubei.crs,drop=True)
ds_clip

以下是数据裁剪前后绘图的对比：

当shp文件含有多个区划时

上面是shp文件仅含有一个地区几何图形的情况，接下来是文件含有多个区划，但是我们仅需要其中某几个区划的数据的情况。同样先读取shp文件：

cn_provinces=gpd.read_file(r'F:\shp文件合集\meteoinfo地图文件\cn_province.shp',encoding='utf-8')

这是一份全国分省级地图大杂烩，假设我们仅需要山东省、安徽省数据，则首先提取山东省、安徽省geometry:

cn_provinces[cn_provinces['NAME'].isin(['山东省','安徽省'])]

然后按照流程开始裁剪数据：

ds2provinces=ds.rio.write_crs('EPSG:4326').rio.clip(cn_provinces[cn_provinces['NAME'].isin(['山东省','安徽省'])].geometry.values,cn_provinces.crs,drop=True)

绘制未裁剪850百帕相对湿度图像：

ax.pcolormesh(ds.longitude.values,ds.latitude.values,ds.r[0,5].values,vmin=0,vmax=120,transform=ccrs.PlateCarree())

绘制已裁剪850百帕相对湿度图像：

ax.pcolormesh(ds2provinces.longitude.values,
              ds2provinces.latitude.values,
              ds2provinces.r[0,5].values,
              vmin=0,vmax=120,transform=ccrs.PlateCarree())

绘制未裁剪850百帕风场：

ax.quiver(ds.longitude.values[::2],ds.latitude.values[::2],ds.u[0,5].values[::2,::2],ds.v[0,5].values[::2,::2],transform=ccrs.PlateCarree(),scale=150)

绘制已裁剪850百帕风场：

ax.quiver(ds2provinces.longitude.values[::2],
          ds2provinces.latitude.values[::2],
          ds2provinces.u[0,5].values[::2,::2],
          ds2provinces.v[0,5].values[::2,::2],
          transform=ccrs.PlateCarree(),scale=150)

小结

总的来说rioxarray由于是基于xarray开发的，所以和xarray的关联性特别好，处理完成的数据也是xarray的数据格式，容易读取。作为xarray官方推荐的取代rasterio的接口的新库包，使用起来是比较方便的。
直接对数据进行筛选，缺点就是绘图的时候会出现狗啃痕迹，痕迹的丑陋程度取决于数据的精度。一个1×1°的数据是肯定比0.25×0.25°的数据更丑的。这种方法适合进行区域计算，而不适合绘图。

为什么强调要写入投影再裁剪

为什么我们要强调裁剪的时候被裁剪的数据需要有投影，而我们的裁剪行政区划也需要投影呢。虽然对于绝大部分公开发行的气象资料来说，栅格数据基本都是PlateCarree的，即数据的经纬度往往是最常见的那种： 但是为了避免可能出现的错误，最好指定投影。同时有些数据确实不是常见的投影模式。例如一份非常规的数字高程tif格式数据。

dem_path=r'F:\BaiduNetdiskDownload\湖北省\1000mhubei.tif'
dem=xr.open_dataset(dem_path)
dem

很明显这不是我们常见的那种经纬度格式，经查，该种投影是阿尔伯斯等面积圆锥投影。虽然绘制的图像确实是湖北，但是经纬度合不上； 这种情况下就只能对原始数据进行投影变换，变成我们常见的投影格式：

dem.rio.reproject('epsg:4326')

绘图恢复正常：