Python+CDO批处理CMIP5模式数据 插值

JOJOweibo

春节后被同学安利了Python，沾上以后深深的被它的简洁所吸引。同时最近需要处理CMIP5模式数据，要对N个模式结果插值到指定网格，Climate Data Operators (CDO)在这方面表现出色，因此，我使用这两种语言处理CMIP5数据，得到指定层，指定时间跨度，差值到指定网格的数据。本人菜鸟一枚，在跌跌撞撞中编好了这个小程序，交卷给气象家园。虽然我编程差，还是希望大家能多交流，当然，也希望能帮到有同样需求的小伙伴。
1）首先描述问题：
组里面的CMIP数据是以model名字为目录存放的（图1），比如ACCESS1-3的所有变量（va,ua,ts...)存放在文件夹/model/ACCESS1-3/中

此外，每个变量切割成了几个不同时间跨度的文件（图2），eg：va变量1850-2005年的数据存放在三个文件中：
va_Amon_ACCESS1-3_historical_r1i1p1_185001-189912.nc
va_Amon_ACCESS1-3_historical_r1i1p1_190001-194912.nc
va_Amon_ACCESS1-3_historical_r1i1p1_200001-200512.nc
另外，每个模式的分辨率不同，我们需要插值到相同的网格上处理。
2）解决思路：
遍历model中的变量（e.g. va), 将上述三个va的文件粘合为一个临时文件（outfiletemp1）。
在outfiletemp1中挑出850hPa，以及需要的时间跨度（eg1900-2005）输出为临时文件(outfiletemp2),
然后对outfiletemp2进行插值，cdo remapbil,输出最终文件outfile就可以了


如下为python代码：

1. # JOJO,IAP,Beijing 2017/09/05 Email:mtjsummer@163.com
   
2. import numpy as np
   
3. from netCDF4 import Dataset
   
4. # #CDO做了一个python 的包，大家pip install cdo 就可以 import了
   
5. from cdo import *
   
6. import os
   
7. cdo = Cdo()
   
8. 
   
9. # there are 41 models data
   
10. models=['bcc-csm1-1', 'CanESM2', 'CNRM-CM5', 'FGOALS-g2', 'FGOALS-s2', \
    
11. 'GFDL-CM3', 'GFDL-ESM2G', 'GFDL-ESM2M', 'GISS-E2-H', 'GISS-E2-R',\
    
12. 'MRI-CGCM3', 'MIROC-ESM', 'MIROC-ESM-CHEM', 'MIROC4h', 'MIROC5',\
    
13.   'NorESM1-M', 'HadCM3', 'HadGEM2-CC', 'HadGEM2-ES', 'inmcm4', \
    
14.   'IPSL-CM5A-LR', 'NorESM1-ME', 'CSIRO-Mk3-6-0', 'MPI-ESM-LR', \
    
15.   'CanCM4', 'ACCESS1-0', 'CCSM4', 'IPSL-CM5A-MR', 'IPSL-CM5B-LR',\
    
16.    'EC-EARTH', 'MPI-ESM-P', 'HadGEM2-AO', 'ACCESS1-3', 'BNU-ESM', \
    
17.    'MPI-ESM-MR', 'CMCC-CM', 'CESM1-FASTCHEM', 'FIO-ESM', 'CESM1-CAM5',\
    
18.     'CESM1-BGC', 'CESM1-WACCM']
    
19. 
    
20. 
    
21. 
    
22. <div class="blockcode"><blockquote># 设置全局变量
    
23. varname='va'
    
24. yrSart=str(1948)
    
25. yrLast=str(2005)
    
26. infolder='/disk_raid5/data/CMIP5/historical/model/'
    
27. outfolder='/disk_raid5/pub/mtj/works/0EAWM_ENSO/CMIP5/'+varname+'/'
    
28. i=0
    
29. 
    
30. for model in models:
    
31. 
    
32.     # select target files---Method 0---------
    
33.     filesdir=os.path.join(infolder, model)
    
34.     print(i,model)
    
35.     os.chdir(filesdir)#改变当前工作目录到指定的路径os.chdir()
    
36.     files=[fname for fname in glob.glob(varname+'_Amon_'+model+'_historical_r1i1p1_*.nc')]
    
37.     files.sort(key=str.lower)#将files按字母顺序排序
    
38.     print(files)
    
39.     # select target files---Method 1---------
    
40.     # filesdir=os.path.join(infloder,models[0])
    
41.     # print(filesdir)
    
42.     # fs2=[d for d in os.listdir(filesdir) if os.path.isfile(os.path.join(filesdir,d)) \
    
43.     # and d.split('_')[0]=='va' and d.split('_')[4]=='r1i1p1']
    
44.     # print(fs2)
    
45.     # 检查model中是否存在varname的数据
    
46.     if len(files)==0:
    
47.         print('there is no ',varname,'data in ',model)
    
48.         continue
    
49.   
    
50.     # 检查outfile是否已经计算过了
    
51.     outfile=outfolder+varname+'_Amon_'+model+'_r1i1p1_'+yrSart+'-'+yrLast+'.nc'
    
52.     if os.path.isfile(outfile):
    
53.             print('the ',outfile,' is already existed')
    
54.             continue
    
55.     # combile files temporarily
    
56.     outfiletemp1=outfolder+varname+'_temp1.nc'
    
57.     cdo.cat(input=files,output=outfiletemp1)
    
58.     print('combile files OK')   
    
59. 
    
60.     # cdo to select level,year (outfiletemp2) and remapbil (outfile)
    
61.     outfiletemp2=outfolder+varname+'_temp2.nc'
    
62.     outfiletemp3=outfolder+varname+'_temp3.nc'
    
63.     cdo.sellevel('85000',input=outfiletemp1,output=outfiletemp2)
    
64.     cdo.selyear('1979/2008',input=outfiletemp2,output=outfiletemp3)
    
65.     cdo.remapbil('r144x73',input =outfiletemp3,output =outfile)
    
66.     print('the monthly',yrSart,'-',yrLast,' r144*73')
    
67.     os.remove(outfiletemp1)
    
68.     os.remove(outfiletemp2)
    
69.     os.remove(outfiletemp2)
    
70.     print('tempfiles removed,',model,'is OK')
    
71.     i=i+1

复制代码

萤火飞侠

请问cdo垂直方向上插值有什么要求吗？我用的ml2pl,92500,50000,2000 ifile ofile，但是显示我No 3D variable with hybrid sigma pressure coordinate found！是什么问题呢

小米虫TcT

♂雨已~○ 发表于 2018-1-21 05:51
楼主请问你是在windows下运行的吗？我在windows下运行总是出现cdo=Cdo（）这出错

请问你解决这个问题了嘛？我也是安装上了cdo但是这句命令总报错

麦芽伙计

♂雨已~○ 发表于 2018-1-21 05:51
楼主请问你是在windows下运行的吗？我在windows下运行总是出现cdo=Cdo（）这出错

请问您的问题解决了吗我也是这个问题出错

JOJOweibo

标题错别字了，插值, 用python的小伙伴，欢迎交流QQ：305897795

贫道敬孔

多谢楼主分享{:5_214:}

小黄-英杰

大神

tiantianyixin

请教一下，我准备用cat函数合并nc数据，错误提示为unsupported file type（library support not complied in）。这是因为没有安装netcdf4扩展包吗？如何安装呢？

明明~~

学到啦~谢楼主

JOJOweibo

tiantianyixin 发表于 2017-5-10 11:50
请教一下，我准备用cat函数合并nc数据，错误提示为unsupported file type（library support not complied i ...

可能是这句有一点点小错误，少了一个空格，cdo无法读入： os.system('cdo select,level=85000,year='+yrSart+'/'+yrLast+outfiletemp1+' '+outfiletemp2)
你先试试改成这样：os.system('cdo select,level=50000,year=1950/2000 '+outfiletemp1+' '+outfiletemp2)
如果有问题再QQ我

JOJOweibo

tiantianyixin 发表于 2017-5-10 11:50
请教一下，我准备用cat函数合并nc数据，错误提示为unsupported file type（library support not complied i ...

推荐使用anaconda安装各种包，

JOJOweibo

tiantianyixin 发表于 2017-5-10 11:50
请教一下，我准备用cat函数合并nc数据，错误提示为unsupported file type（library support not complied i ...

推荐使用anaconda安装各种包，

felin

厉害了，才学Python就这么666