ERA5温度、风速风向处理

pro READ_ERA5_DAILY_wind
  ;该程序用来处理ERA5的skin temperature数据
  file_ID = NCDF_OPEN('E:\heatflux\ERA5\ERA5_20190420_1600.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
  ;'G:\Lead\data\validation\20180222T200801\ERA5\adaptor.mars.internal-1694329666.815577-13956-2-95bf6c78-4b8c-4620-bac4-689b9456e475.nc'
  outfilepath='E:\heatflux\ERA5\'
  ;'G:\Lead\data\validation\20180222T200801\ERA5\'
  ;'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,'u10')   ;  'msl'           ;g0_lon_3
  NCDF_VARGET, file_ID, mslid, u10          ;获取 land维变量 int
  nmsl = N_ELEMENTS(u10)
  
  vlid = NCDF_VARID(file_ID,'v10')   ;  'msl'           ;g0_lon_3
  NCDF_VARGET, file_ID, vlid, v10          ;获取 land维变量 int
  
  d2mid = NCDF_VARID(file_ID,'d2m')   ;  'msl'           ;g0_lon_3
  NCDF_VARGET, file_ID, d2mid, d2m          ;获取 land维变量 int
  
  t2mid = NCDF_VARID(file_ID,'t2m')   ;  'msl'           ;g0_lon_3
  NCDF_VARGET, file_ID, t2mid, t2m          ;获取 land维变量 int
;  nmsl = N_ELEMENTS(msl)

  ncdf_attget,file_ID,mslid,'scale_factor',ua
  ncdf_attget,file_ID,mslid,'add_offset',ub
  
  ncdf_attget,file_ID,vlid,'scale_factor',va
  ncdf_attget,file_ID,vlid,'add_offset',vb
  
  ncdf_attget,file_ID,d2mid,'scale_factor',da
  ncdf_attget,file_ID,d2mid,'add_offset',db

  ncdf_attget,file_ID,t2mid,'scale_factor',ta
  ncdf_attget,file_ID,t2mid,'add_offset',tb

;  print,max(msl),min(msl),mean(msl)
  if size(latitude,/N_DIMENSIONS) eq 1 then begin ;如果经纬度为一维的，就要建立成二维的
    templa=fltarr(nlon,nlat)
    templo=fltarr(nlon,nlat)
    ;      templa=MAKE_ARRAY(nlat,nlon,/FLOAT,VALUE =latitude[*])
    for i=0,nlat-1 do begin
      templa[*,i]=latitude[i]
    endfor
    for i=0,nlon-1 do begin
      templo[i,*]=longitude[i]
    endfor

    latitude=templa
    longitude=templo
  endif

  WRITE_TIFF,outfilepath+'ERA5_lon.tif',longitude,/float
  WRITE_TIFF,outfilepath+'ERA5_lat.tif',latitude,/float
  deg=180.0/!DPI
  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

      u10=u10[*,*]*1.0*ua+ub;*double(a)+b
      v10=v10[*,*]*1.0*va+vb
      outtmp=180.0 + ATAN(u10, v10)*deg ;角度
;      注意，这里计算风向角度用的公式是按照气象学上的定义。
;      气象上定义正北方向为0°（即风从北吹向南，是数学坐标系中的 -90°）， 顺时针转动角度增大。
;      风向Dir=0°（或360°）, u=0, v<0, 正北风；
;      风向Dir=90°, u<0, v=0，正东风；
;      风向Dir=180°, u=0, v>0，正南风；
;      风向Dir=270°, u>0, v=0，正西风。
      d2m=d2m*da+db
      t2m=t2m*ta+tb

      windspeed=SQRT(u10^2*v10^2) ;风速Wind (m/s)
;      index=where(u10 eq -32767S,count)
;;      outtmp[*,*,*]=outtmp[*,*,*]*a+b
;      if count gt 0 then begin
;        outtmp[index]=-1 ;无效值被复位为-1
;      endif

      
      WRITE_TIFF,outfilepath+'ERA5_winddirection_'+strtrim(String(year),2)+'_'+strtrim(String(month),2)+'_'+strtrim(String(day),2)+'_'+strtrim(String(hour),2)+'.tif',outtmp,/float; , geotiff=geoinfo;存储数据,会自动覆盖同名数据
      WRITE_TIFF,outfilepath+'ERA5_windspeed_'+strtrim(String(year),2)+'_'+strtrim(String(month),2)+'_'+strtrim(String(day),2)+'_'+strtrim(String(hour),2)+'.tif',windspeed,/float; , , PLANARCONFIG=2geotiff=geoinfo;
      WRITE_TIFF,outfilepath+'ERA5_d2m_'+strtrim(String(year),2)+'_'+strtrim(String(month),2)+'_'+strtrim(String(day),2)+'_'+strtrim(String(hour),2)+'.tif',d2m,/float; , geotiff=geoinfo;存储数据,会自动覆盖同名数据
      WRITE_TIFF,outfilepath+'ERA5_t2m_'+strtrim(String(year),2)+'_'+strtrim(String(month),2)+'_'+strtrim(String(day),2)+'_'+strtrim(String(hour),2)+'.tif',t2m,/float; , , PLANARCONFIG=2geotiff=geoinfo;

  endfor


end

后面增加个地理校正批处理