赞
踩
首先是我的数据下载方面,从官网上直接下载没有使用代码,因此所有下载的数据都放到一个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
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。