在对全球环境进行研究中，我们总会遇到欧洲中期预报中（European Centre for Medium-Range Weather Forecasts）等发布的气象再分析数据，例如，我在ESA官网上下的全球逐日开放水体蒸散数据，都是以GRIB格式结尾的。(以ubuntu 18.04为例)

首先要下载

https://confluence.ecmwf.int//display/ECC/Releases

https://confluence.ecmwf.int/download/attachments/45757960/eccodes-2.18.0-Source.tar.gz?api=v2 （下载最新版文件）

mkdir build

tar -xvf eccodes-2.18.0-Source.tar.gz

cd build

cmake ../eccodes-2.18.0-Source -DCMAKE_INSTALL_PREFIX=/usr/local -DCMAKE_C_COMPILER=gcc -DCMAKE_Fortran_COMPILER=gfortran-9

ldconfig （这个命令是必须的，不然找不到libecodes.so错误就会出现，前面那个 gfortan-9记得加上去，因为gfortan版本也在变，用哪个版本的fortan就加哪个数字在后面）

下载cdo最新版本

https://code.mpimet.mpg.de/projects/cdo

wget https://code.mpimet.mpg.de/attachments/22390/cdo-1.9.9rc6.tar.gz

tar -xvf cdo-1.9.9rc6.tar.gz

cd cdo-1.9.9rc6

mkdir build

./configure --with-netcdf --with-eccodes=/usr/local

make && make install

然后可以使用cdo info **.grib查看文件信息是否可以读取

在python下面使用可以使用参考

https://github.com/ecmwf/cfgrib 使用xarray来读grib文件，但有个问题是同样速度会很慢的，这里介绍另外一种方式。因为，我们研究区域通常会很小，可以考虑使用

cdo -sellonlatbox 命令 （不多说了，直接上代码）

import geopandas as gpd #使用geopandas读取全球湖泊水库矢量

df = gpd.read_file(r"./shapefilegallery/Export_Output.shp")

df.crs = {'init': 'epsg:4326', 'no_defs': True}

df = df.to_crs("EPSG:4326") # 都转换成EPSG:4326

lake_reserviro = utility.getmaxreserviors_lakes()

#print(lake_reserviro[1])

lakes = lake_reserviro[1][:] # 假定这个是所有湖泊编号

input_grib_file = **.grib #假定这个是要读取的XX.grib文件

for specify_which_reservior in lakes:

the_geometry = df.loc[df['GLWD_ID']==specify_which_reservior,'geometry'].iloc[0]

print(the_geometry.bounds)

lowlon = the_geometry.bounds[0]

uplong = the_geometry.bounds[2]

lowlat = the_geometry.bounds[1]

uplat = the_geometry.bounds[3]

output_name = os.path.join("./",str(specify_which_reservior)+"xx"+ ".grib")

command = " cdo -sellonlatbox," + str(lowlon) + ","+ str(uplong) + ","+str(lowlat) + ","+str(uplat)

command = command + " " + input_grib_file+" "+output_name

os.system(command) #用该命令输出指定湖泊范围内的grid文件

接着，我们可以使用grib_get_data 获取指定时间段的气象数据，然后将结果存放到csv文件中

year = the_time.year

month = the_time.month

day = the_time.day

the_query_date = str(year)+str(month).rjust(2,'0')+str(day).rjust(2,'0')

output_name = os.path.join("./",the_query_date+"_precipitation"+ ".csv")

command = " grib_get_data " + " -w dataDate=" + the_query_date + " " + the_grib_file + " > "

command = command + output_name

os.system(command)

这里的the_query_date 一般采用年月日形式传递，例如20090101

后面，我们会继续分享，如何得到指定湖泊区域100公里范围内的气象再分析数据。本文来自于善睐科技-黄伟。谢谢大家持续关注我们~~~~