IDL 处理ERA5的NC格式的温度及海平面气压数据

亦凉 2022-12-25 13:51 233阅读 0赞

首先是我的数据下载方面，从官网上直接下载没有使用代码，因此所有下载的数据都放到一个NC文件中了（当然2m温度和气压是分开的）

我主要下载的区域是北纬30-90，每日只下载四个时间，0 6 12 18，也只下载6个月的，1-3、10-12月

nc文件可以尝试使用panoly打开看看，每个月和每个时间是按顺序排好的

我利用以下的代码来处理2m温度和海平面气压数据

值得注意的点是，在读取数据时，还要读取校正系数和偏差，即下面的a和b变量，并且注意读取时数据格式是int，运算时需要先对数组乘以1.0变成float类型的，不然运算就按整数进行，会出错。

里面地理坐标的设置，我是通过先将nc数据放入arcmap里面，直接导出一个tif，最后将这个tif的地理信息赋给所有文件。如果有直接通过IDL来获取的方法，欢迎留言

pro ERA5_Daily_SLP_read
    ;该程序用来处理2015-2019年ERA50的海平面气压数据
        file_ID = NCDF_OPEN('H:\IceClassification\sentinel\Lancaster\ERA5\sea surface temperature\30_90_2mt_2015_2019_4time_1_3_10_12month.nc',/NOWRITE) ;'H:\IceClassification\sentinel\Lancaster\ERA5\sea level pressure\30_90_SLP_2015_2019_4time_1_3_10_12month.nc'  ;open netCDF file for READ only
        outfilepath='H:\IceClassification\sentinel\Lancaster\ERA5\sea surface temperature\2mt_tif\';'H:\IceClassification\sentinel\Lancaster\ERA5\sea level pressure\sea_level_pressure_tif\';;
        geo=read_tiff('H:\IceClassification\mask\LS\ERA5_30_geo.tif',geotiff=geoinfo);
        Tag = NCDF_INQUIRE(file_ID)
        timeid = NCDF_VARID(file_ID,'time')             ;initial_time0initial_time0_hours  initial_time0
        NCDF_VARGET, file_ID, timeid, time              ;获取 时间维变量 long
        nt = N_ELEMENTS(time)
        ;print,time
        latid = NCDF_VARID(file_ID,'latitude')               ;g0_lat_2 flo
        NCDF_VARGET, file_ID, latid, latitude           ;获取 时间维变量
        nlat = N_ELEMENTS(latitude)                     ;纬向维长度 float
        ;print,'dsklfjdsjfkldjsfkljl'+latitude
        lonid = NCDF_VARID(file_ID,'longitude')                ;g0_lon_3
        NCDF_VARGET, file_ID, lonid, longitude           ;获取 时间维变量
        nlon = N_ELEMENTS(longitude)                     ;径向维长度 float
      
        mslid = NCDF_VARID(file_ID,'t2m')   ;  'msl'           ;g0_lon_3
        NCDF_VARGET, file_ID, mslid, msl          ;获取 land维变量 int
        nmsl = N_ELEMENTS(msl)
      
        ncdf_attget,file_ID,mslid,'scale_factor',a
        ncdf_attget,file_ID,mslid,'add_offset',b
    ;    CALDAT,(time[nt-1]+julday(01,01,1900,0,0,0)*24.0)/24.0,month,day,year,hour,minu;,sec
    ;    print,month,day,year,hour,minu,sec
      print,max(msl),min(msl),mean(msl)
        for t=0,nt-1 do begin
          CALDAT,(time[t]+julday(01,01,1900,0,0,0)*24.0)/24.0,month,day,year,hour,minu;,sec
           print,month,day,year,hour,minu;,sec
           if hour eq 0 then begin
              outtmp=msl[*,*,t:t+3]*1.0;*double(a)+b
              index=where(outtmp eq -32767S,count)
              outtmp[*,*,*]=outtmp[*,*,*]*a+b
              if count gt 0 then begin
                outtmp[index]=0 ;无效值背负为0
              endif
              WRITE_TIFF,outfilepath+'ERA50_2mt_'+strtrim(String(year),2)+'_'+strtrim(String(month),2)+'_'+strtrim(String(day),2)+'.tif',outtmp, PLANARCONFIG=2,/float, geotiff=geoinfo; ;存储数据,会自动覆盖同名数据
           endif
           
        endfor
        
    
    end