下载了一些模型包括以下:
AWI-CM-1-1-MR |
ACCESS-CM2 |
BCC-CSM2-MR |
CESM2 |
CESM2-WACCM |
CMCC-CM2-SR5 |
CMCC-ESM2 |
CNRM-CM6-1 |
CNRM-ESM2-1 |
CanESM5 |
EC-Earth3 |
EC-Earth3-Veg |
EC-Earth3-Veg-LR |
IPSL-CM6A-LR |
MIROC-ES2L |
MIROC6 |
MPI-ESM1-2-HR |
MPI-ESM1-2-LR |
MRI-ESM2-0 |
NorESM2-LM |
NorESM2-MM |
EC-Earth3-CC |
KIOST-ESM |
我自己下载的是siconc和sithick两个变量,但是本程序已经可以根据文件名自动判别变量名称了,而且存在SIMon,和SIday两种frequency的。
所以要求程序能够自动处理不同模型、变量不同、频率不同的情况。
程序只把数据提取成envi标准格式,并且也存了envi的hdr,文件名需要保留除日期之外的全部,再加上日期/月份,存储的文件是月/日的。
最后也保存了lon和lat
我觉得比较难做的主要是,各种日期的识别,因为有考虑闰年的,不考虑闰年的之类,因此都单独用了函数来计算。
pro read_CMIP6_model_data
infilepath = 'H:\mission\LNG\Model\MPI-ESM1-2-LR\ssp585\daily\sicon'
outfilepath = 'H:\mission\LNG\Model\MPI-ESM1-2-LR\ssp585\daily\sicon\hdr\'
nv = envi(/HEADLESS)
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)
NCDF_ATTGET,file_ID, timeid,'calendar',timeinfo
timeinfo=string(timeinfo)
NCDF_ATTGET,file_ID, timeid,'units',timeunits
timeunits=string(timeunits)
timeunits=STRSPLIT(timeunits, ' ', /EXTRACT)
;print,time
latid = NCDF_VARID(file_ID,'lat') ;g0_lat_2
if latid eq -1 then begin
latid = NCDF_VARID(file_ID,'latitude')
if latid eq -1 then begin
latid = NCDF_VARID(file_ID,'nav_lat')
endif
endif
NCDF_VARGET, file_ID, latid, latitude ;获取 维度变量
nlat = N_ELEMENTS(latitude) ;纬向维长度
;print,'dsklfjdsjfkldjsfkljl'+latitude
lonid = NCDF_VARID(file_ID,'lon')
if lonid eq -1 then begin
lonid = NCDF_VARID(file_ID,'longitude')
if lonid eq -1 then begin
lonid = NCDF_VARID(file_ID,'nav_lon')
endif
endif ;g0_lon_3
NCDF_VARGET, file_ID, lonid, longitude ;获取 经度变量
nlon = N_ELEMENTS(longitude) ;径向维长度
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
var_ID = NCDF_VARID(file_ID,str[0] ) ;获取变量IDvar_result.Name siconc'sithick'
NCDF_VARGET, file_ID, Var_ID, var_Data ;获取变量数据
NCDF_ATTGET,file_ID, var_ID,'missing_value',varinfo
varinfo=float(varinfo)
timeunits=STRSPLIT(timeunits[2], '-', /EXTRACT)
year=long(timeunits[0]);1850;long(STRMID(str[6],0,4))
month=long(timeunits[1]);long(STRMID(str[6],4,2))
day=long(timeunits[2])
NCDF_CLOSE, file_ID
print,' All parameters of file ' + thesefiles[fidx] + ' have been writen out!'
index=where(latitude[*,*] ge 47,count) ;我只要纬度>47的数据
ind = ARRAY_INDICES(latitude,index)
; 输出经纬度
dim=size(latitude[index],/DIMENSIONS)
if N_elements(dim) eq 1 then begin
dim=[dim,1]
endif
envi_setup_head,FNAME=infilepath+'\lat',NS=dim[0] ,NL=dim[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=dim[0] ,NL=dim[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
;通过时间变量的计算判定日期
for x=0,nt-1 do begin
if timeinfo eq '360_day' then begin
tempresult=calendar_360(year,month,day,time[x])
yr=tempresult[0]
mon=tempresult[1]
dy=tempresult[2]
endif else if timeinfo eq 'noleap' then begin
tempresult=calendar_noleap(year,month,day,time[x])
yr=tempresult[0]
mon=tempresult[1]
dy=tempresult[2]
endif else begin
caldat,julday(month,day,year)+time[x]-0.5,mon,dy,yr;这个是proleptic_gregorian或者儒略历(Julian) thick需要在filename后面添加这个
; print,mon,dy,yr
; gregorian and julday For dates between 1 Jan CE and 15 Oct 1582 the two calendar systems differ by up to 10 days. For dates on or after 15 Oct 1582 the two calendar systems are identical.
endelse
if yr lt 1980 then begin
continue
endif
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)')
if mon eq 2 and dy eq 29 then begin ;不要闰年的2月29号数据
continue
endif
if yr gt 2059 then begin
continue
endif
endif
; temp=[]
if size(var_Data,/N_DIMENSIONS) eq 3 then begin ;这是由于变量不是像AWI,而是两维+1维时间
dim=size(var_Data,/DIMENSIONS)
temp=fltarr(dim[0],dim[1])
temp[*,*]=var_Data[*,*,x]
temp=temp[index];var_Data[ind[1],ind[0],x]
endif else begin
temp=var_Data[index,x]
endelse
miss=where(temp eq varinfo,count)
if count gt 0 then begin
temp[miss]=-1
endif
; temp1=var_Data[ind[0],ind[1],x]
; print,max(temp),mean(temp),var_Data[1234,137,x];,temp[100,37]
dim=size(temp,/DIMENSIONS)
if N_elements(dim) eq 1 then begin
dim=[dim,1]
endif
envi_setup_head,FNAME=outfilename,NS=dim[0] ,NL=dim[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
;消除内存
temp=!Null
delvar,temp
print,outfilename
endfor
; 输出经纬度
dim=size(latitude[index],/DIMENSIONS)
if N_elements(dim) eq 1 then begin
dim=[dim,1]
endif
envi_setup_head,FNAME=infilepath+'\lat',NS=dim[0] ,NL=dim[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=dim[0] ,NL=dim[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
var_Data=!Null
delvar,var_Data
print,thesefiles[fidx],' is finished'
ENDFOR
nv.close
end
function calendar_360,yr,mon,dy,time
tdy=dy+time mod 30
tmon=long(time/30)
tyr=yr+long(tmon/12)
tmon=mon+tmon mod 12
return,[tyr,tmon,tdy]
end
FUNCTION calendar_noleap,yr,mon,dy,time
;不考虑闰年
;
months_days = [31,28,31,30,31,30,31,31,30,31,30,31]
time=time+0.5
; IF ISLEAPYEAR(year) EQ 1 THEN months_days[1]=28;29 ;注意我这边由于闰年的28都被排除了,所以DOY不计算29的情况
tyr=yr+long(time/365)
tdy=time mod 365
for month=0,12-1 do begin
if tdy le total(months_days[0:month]) then begin
break
endif
endfor
tmon=month+1
if month gt 0 then begin
tdy=tdy-total(months_days[0:month-1])
endif
RETURN, [tyr,tmon,tdy]
END