下载了一些模型包括以下:
AWI-CM-1-1-MR |
CMCC.CMCC-CM2-HR4 |
CNRM-CERFACS.CNRM-CM6-1-HR |
NOAA-GFDL.GFDL-CM4 |
MOHC.HadGEM3-GC31-MM |
下载的是siconc和sithick两个变量,而且存在SIMon,和SIday两种frequency的。
所以要求程序能够自动处理不同模型、变量不同、频率不同的情况。
程序只把数据提取成envi标准格式,并且也存了envi的hdr,文件名需要保留除日期之外的全部,再加上日期/月份,存储的文件是月/日的。
最后也保存了lon和lat
pro read_CMIP6_model_data
infilepath = 'H:\Model\AWI-CM-1-1-MR\historical\scon\daily'
outfilepath = 'H:\Model\AWI-CM-1-1-MR\historical\scon\daily\hdr\'
nv = envi(/HEADLESS)
; geo=read_tiff('H:\mission\LNG\mask\smos_geo2.tif',geotiff=geoinfo);
; mask=read_tiff('H:\mission\LNG\mask\smos_landmask.tif');
CD,infilepath
thesefiles = FILE_SEARCH('*.nc',count = numFiles)
print,'number of files found:',numFiles
;处理数据
FOR fidx=0, numFiles-1 DO BEGIN; ;不要忘记修改这个
filename=file_basename(thesefiles[fidx],'.nc')
str=STRSPLIT(filename, '_', /EXTRACT) ;siconc_SIday_AWI-CM-1-1-MR_historical_r1i1p1f1_gn_18500101-18501231
;获取的是一期的数据
file_ID = NCDF_OPEN(theseFiles[fidx],/NOWRITE) ;open netCDF file for READ only
Tag = NCDF_INQUIRE(file_ID)
timeid = NCDF_VARID(file_ID,'time') ;initial_time0initial_time0_hours initial_time0
NCDF_VARGET, file_ID, timeid, time ;获取 时间维变量
nt = N_ELEMENTS(time)
;print,time
latid = NCDF_VARID(file_ID,'lat') ;g0_lat_2
if latid eq -1 then begin
latid = NCDF_VARID(file_ID,'latitude')
endif
NCDF_VARGET, file_ID, latid, latitude ;获取 时间维变量
nlat = N_ELEMENTS(latitude) ;纬向维长度
;print,'dsklfjdsjfkldjsfkljl'+latitude
lonid = NCDF_VARID(file_ID,'lon')
if latid eq -1 then begin
latid = NCDF_VARID(file_ID,'longitude')
endif ;g0_lon_3
NCDF_VARGET, file_ID, lonid, longitude ;获取 时间维变量
nlon = N_ELEMENTS(longitude) ;径向维长度
var_ID = NCDF_VARID(file_ID,str[0] ) ;获取变量IDvar_result.Name siconc'sithick'
NCDF_VARGET, file_ID, Var_ID, var_Data ;获取变量数据
;大于等于0为海冰有效值,缺失值为-0.2,陆地为-0.1
year=long(STRMID(str[6],0,4))
month=long(STRMID(str[6],4,2))
day=1
if str[1] eq 'SIday' then begin
day=long(STRMID(str[6],6,2))
endif
NCDF_CLOSE, file_ID
print,' All parameters of file ' + thesefiles[fidx] + ' have been writen out!'
index=where(latitude ge 47,count) ;只要纬度>47的数据
for x=0,nt-1 do begin
caldat,julday(month,day,year)+time[x]-0.5,mon,dy,yr;thick需要在filename后面添加这个
outfilename = outfilepath +String(yr,FORMAT='(I04)')+'\'+strjoin(str[0:5],'_')+'_'+String(yr,FORMAT='(I04)')+String(mon,FORMAT='(I02)');+String(dy,FORMAT='(I02)');STRCOMPRESS(outfilepath + 'smos_sea_ice_thickness_'+year+'_'+month+'_'+day+'.tif',/remove_all)
;存储数据
if str[1] eq 'SIday' then begin
outfilename = outfilename+String(dy,FORMAT='(I02)')
endif
temp=var_Data[index,x]
envi_setup_head,FNAME=outfilename,NS=count,NL=1,NB=1,DATA_TYPE=size(temp,/type),INTERLEAVE=0,/OPEN,/WRITE ;4为float
openw,lunou,outfilename,/get_lun
writeu,lunou,temp;
close,lunou
free_lun,lunou
;消除内存
delvar,temp
print,outfilename
endfor
; 输出经纬度
envi_setup_head,FNAME=infilepath+'\lat',NS=count,NL=1,NB=1,DATA_TYPE=size(latitude,/type),INTERLEAVE=0,/OPEN,/WRITE ;5为double
openw,lunou,infilepath+'\lat',/get_lun
writeu,lunou,latitude[index]
close,lunou
free_lun,lunou
envi_setup_head,FNAME=infilepath+'\lon',NS=count,NL=1,NB=1,DATA_TYPE=size(longitude,/type),INTERLEAVE=0,/OPEN,/WRITE
openw,lunou,infilepath+'\lon',/get_lun
writeu,lunou,longitude[index]
close,lunou
free_lun,lunou
; var_Data[where(~finite(var_Data))] = -1
print,thesefiles[fidx],' is finished'
ENDFOR
nv.close
end