CMIP6 数据——批量nc转tif-温度tas（6）
pro CMIP6_batch_tas
netcdf_path = 'G:\CMIP6_\'
;;;model_name
;
;model_name = 'ACCESS-ESM1-5'
;var_name = ['tos','longitude','latitude']
;
;model_name = 'CESM2-WACCM'
;var_name = ['tos','lon','lat']
;
;model_name = 'GFDL-CM4'
;var_name = ['tos','lon','lat']
;
;model_name = 'GFDL-ESM4'
;var_name = ['tos','lon','lat']
;
;model_name = 'IPSL-CM6A-LR'
;var_name = ['tos','nav_lon','nav_lat']
;
;model_name = 'MIROC6'
;var_name = ['tos','longitude','latitude']
;
;model_name = 'MPI-ESM1-2-HR'
;var_name = ['tos','longitude','latitude']
;
;model_name = 'MRI-ESM2-0'
;var_name = ['tos','longitude','latitude']
;model_name = 'NorESM2-LM'
;var_name = ['tas','lon','lat']
;model_name = 'BCC-CSM2-MR'
;var_name = ['tas','lon','lat']
;
;model_name = 'CMCC-ESM2'
;var_name = ['tas','lon','lat']
;
;model_name = 'ACCESS-CM2'
;var_name = ['tas','lon','lat']
;model_name = 'CanESM5'
;var_name = ['tas','lon','lat']
;model_name = 'CESM2-WACCM'
;var_name = ['tas','lon','lat']
;model_name = 'CMCC-CM2-SR5'
;var_name = ['tas','lon','lat']
;model_name = 'FGOALS-g3'
;var_name = ['tas','lon','lat']
;model_name = 'MIROC6'
;var_name = ['tas','lon','lat']
;model_name = 'MPI-ESM1-2-HR'
;var_name = ['tas','lon','lat']
;model_name = 'MPI-ESM1-2-LR'
;var_name = ['tas','lon','lat']
;
model_name = 'MRI-ESM2-0'
var_name = ['tas','lon','lat']
file_arr = file_search(netcdf_path,'tas*'+model_name+'*.nc')
nfile = N_elements(file_arr)
for fi = 0,nfile-1 do begin
file_path = file_arr[fi]
file_name = FILE_BASENAME(file_path)
str_split = strsplit(file_name,'_',/EXTRACT)
perfname = str_split[2]+'_'+str_split[3]+'_'+'tas_'
start_date = strmid(str_split[6],0,6)
str_split = strsplit(file_name,'.',/EXTRACT)
data_name = str_split[0]+'.dat'
geo_data_name = 'geo_'+str_split[0]+'.dat'
writepath = 'G:\CMIP6_tas\GEOTIFF\'+ model_name+'\'
glt_path = writepath + 'glt\'
a = read_ncfile_to_envigeotiff(file_path,glt_path,data_name,geo_data_name,var_name)
b = read_envigeotiff_to_geotiff(model_name,glt_path+geo_data_name,writepath,start_date,perfname)
temp=0;
print,fi,data_name
endfor
end
function read_ncfile_to_envigeotiff,file_path,glt_path,data_name,geo_data_name,var_name
; e = ENVI()
COMPILE_OPT idl2
compile_opt strictarr
ENVI,/RESTORE_BASE_SAVE_FILES
ENVI_BATCH_INIT ,LOG_FILE = 'batch.log'
e=ENVI(/headless);处理ENVI命令，但不显示envi窗体
; envi_batch_status_window,/on ;/on参数：显示处理进度条
if file_test(glt_path, /directory) eq 0 then file_mkdir,glt_path;输出路径没有即新建立
writepath = glt_path;
filepath = file_path
nc_id = NCDF_OPEN(filepath,/nowrite); 1.读取文件的id
var_id = NCDF_VARID(nc_id,var_name[0]); 2.读取变量的id
NCDF_VARGET,nc_id,var_id,tos; 3.读取变量值
var_id = NCDF_VARID(nc_id,var_name[1]); 2.读取变量的id
NCDF_VARGET,nc_id,var_id,lon_temp; 3.读取变量值
var_id = NCDF_VARID(nc_id,var_name[2]); 2.读取变量的id
NCDF_VARGET,nc_id,var_id,lat_temp; 3.读取变量值
NCDF_CLOSE,nc_id; 关闭
temp= 0 ;
nt=N_elements(tos[0,0,*]);获取时间数
ns=N_elements(tos[*,0,0]);获取列数
nl=N_elements(tos[0,*,0]);获取行数
longitude = fltarr(ns,nl)
latitude = fltarr(ns,nl)
for i=0,ns-1 do begin
latitude[i,*] = lat_temp[*]
endfor
for j=0,nl-1 do begin
longitude[*,j] = lon_temp[*]
endfor
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;存储经纬度数据并构建glt查找表
if file_test(writepath+'RefLon.dat') eq 0 and file_test(writepath+'RefLat.dat') eq 0 then begin
RefLon = ENVIRaster(longitude,uri=writepath+'RefLon.dat')
RefLat = ENVIRaster(latitude ,uri=writepath+'RefLat.dat')
RefLon.save
RefLat.save
RefLon_id = ENVIRasterToFID(RefLon)
RefLat_id = ENVIRasterToFID(RefLat)
endif
if file_test(writepath+'RefLon.dat') eq 1 and file_test(writepath+'RefLat.dat') eq 1 then begin
ENVI_OPEN_DATA_FILE,writepath+'RefLon.dat',r_fid=RefLon_id
envi_file_query,RefLon_id,dims=dims_x,nb=nb_lon
x_pos=lindgen(nb_lon)
ENVI_OPEN_DATA_FILE,writepath+'RefLat.dat',r_fid=RefLat_id
envi_file_query,RefLat_id,dims=dims_y,nb=nb_lat
y_pos=lindgen(nb_lat)
endif
i_proj=ENVI_PROJ_CREATE(/geographic, datum="WGS-84")
o_proj=ENVI_PROJ_CREATE(/geographic, datum="WGS-84")
if file_test(writepath+'glt.tif') eq 0 then begin
ENVI_DOIT,'ENVI_GLT_DOIT', rotation=0, pixel_size = 1, DIMS = dims,I_PROJ= i_proj, O_PROJ=o_proj, out_name = writepath+'glt.tif', $
R_FID=glt_fid, X_FID=RefLon_id, X_POS=x_pos[0], Y_FID=RefLat_id, Y_POS=y_pos[0]
endif
if file_test(writepath+'glt.tif') eq 1 then begin
ENVI_OPEN_DATA_FILE,writepath+'glt.tif',r_fid=glt_fid
ENVI_FILE_QUERY, glt_fid, dims=dims
endif
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;存储tos数据并做glt校正
tos_name = writepath + data_name
geotos_name = writepath + geo_data_name
if file_test(tos_name) eq 0 then begin
Refdata = ENVIRaster(tos,uri=tos_name)
Refdata.save
Refdata_id = ENVIRasterToFID(Refdata)
endif
if file_test(tos_name) eq 1 then begin
ENVI_OPEN_DATA_FILE,tos_name,r_fid=Refdata_id
envi_file_query,Refdata_id,dims=dims_data,nb=nb_data
endif
ENVI_DOIT,'ENVI_GEOREF_FROM_GLT_DOIT', BACKGROUND=0,FID=Refdata_id,GLT_FID=glt_fid,OUT_NAME=geotos_name,POS=[0:nb_data-1]
;ENVI_DOIT, 'ENVI_GEOREF_FROM_GLT_DOIT', BACKGROUND=-999, FID=Refdata_id,GLT_DIMS=glt_dims, pos=pos, GLT_FID=glt_fid,$
; out_name=geotos_name, r_fid=out_fid
temp=0;
;;;;;;释放内存
fids = envi_get_file_ids();获取当前内存中的所有文件的fid
size = size(fids);获取数组的大小
length = size[1]
;循环释放内存中的文件
for i = 0L, length-1 do begin
envi_file_mng,id = fids[i],/remove
endfor
end
function read_envigeotiff_to_geotiff,model_name,readpath,writepath,start_date,perfname
;;;;;所有的数据统一为-180到180，90到-90，分辨率为1度
COMPILE_OPT idl2
compile_opt strictarr
ENVI,/RESTORE_BASE_SAVE_FILES
ENVI_BATCH_INIT ,LOG_FILE = 'batch.log'
e=ENVI(/headless);处理ENVI命令，但不显示envi窗体
start_year = double(strmid(start_date,0,4))
geotiff = {$
MODELPIXELSCALETAG: [1.0000000000000000, 1.0000000000000000, 0.00000000000000000],$
MODELTIEPOINTTAG: [0.00000000000000000, 0.00000000000000000, 0.00000000000000000, -180, 90, 0.00000000000000000],$
GTMODELTYPEGEOKEY: 2,$
GTRASTERTYPEGEOKEY: 1,$
GEOGRAPHICTYPEGEOKEY: 4326,$
GEOGCITATIONGEOKEY: "GCS_WGS_1984",$
GEOGANGULARUNITSGEOKEY: 9102,$
GEOGSEMIMAJORAXISGEOKEY: 6378137.0000000000,$
GEOGINVFLATTENINGGEOKEY: 298.25722356300003}
ENVI_OPEN_DATA_FILE,readpath,r_fid=data_id
envi_file_query,data_id,dims=dims_data,nb=nb_data
o_proj = ENVI_GET_PROJECTION(FID =data_id)
raster = e.OpenRaster(readpath)
upleft = raster.SPATIALREF.tie_point_map
out_arr = fltarr(nb_data,360,180)
;依次输出月份
writepath1 = writepath+'monthly\'
if file_test(writepath1, /directory) eq 0 then file_mkdir,writepath1;输出路径没有即新建立
for fi = 0,nb_data-1 do begin
;print,fi
data = envi_get_data(fid = data_id, dims = dims_data, pos = fi)
month = fi mod 12 +1.0
; if month eq 1 or month eq 3 or month eq 5 or month eq 7 or month eq 8 or month eq 10 or month eq 12 then begin
; data = data*60.0*60.0*24*31.0
; end
; if month eq 4 or month eq 6 or month eq 9 or month eq 11 then begin
; data = data*60.0*60.0*24*30.0
; end
; if month eq 2 then begin
; data = data*60.0*60.0*24*29.0
; end
ns=N_elements( data[*,0]);获取列数
nl=N_elements( data[0,*]);获取行数
for i = 0,ns-1 do begin
for j = 0,nl-1 do begin
if (model_name eq 'MRI-ESM2-0') then begin
if j gt 0 and j lt nl-1 and data[i,j] eq 0 then data[i,j] = (data[i,j-1] + data[i,j+1])/2.0
if j eq 0 and data[i,j] eq 0 then data[i,j] = data[i,j+1]
if i eq nl-1 and data[i,j] eq 0 then data[i,j] = data[i,j-1]
endif
if data[i,j] lt -1000 then data[i,j] = -999
endfor
endfor
data = data - 273.15
case model_name of
'ACCESS-ESM1-5': begin
out_arr[fi, 0:179,0:nl-1] = data[180:359,0:nl-1]
out_arr[fi,180:359,0:nl-1] = data[0:179,0:nl-1]
end
; 'CESM2-WACCM': begin
; out_arr[fi, 0:179,9:nl-1+9] = data[179:358,0:nl-1]
; out_arr[fi,181:359,9:nl-1+9] = data[0:178,0:nl-1]
; out_arr[fi,180,9:nl-1+9] = data[0,0:nl-1]
; end
'GFDL-CM4': begin
out_arr[fi, 0:239,0:nl-1] = data[120:359,0:nl-1]
out_arr[fi,240:359,0:nl-1] = data[0:119,0:nl-1]
end
'GFDL-ESM4': begin
out_arr[fi, 0:239,0:nl-1] = data[120:359,0:nl-1]
out_arr[fi,240:359,0:nl-1] = data[0:119,0:nl-1]
end
'IPSL-CM6A-LR': begin
out_arr[fi, 1:359,0:nl-1] = data[0:358,0:nl-1]
out_arr[fi, 0,0:nl-1] = data[ 0,0:nl-1]
out_arr[fi, 253,0:nl-1] = data[ 251,0:nl-1]
end
; 'MIROC6': begin
; out_arr[fi, 0:179,0:nl-1] = data[180:359,0:nl-1]
; out_arr[fi,180:359,0:nl-1] = data[0:179,0:nl-1]
; end
; 'MPI-ESM1-2-HR': begin
; out_arr[fi, 0:179,1:nl] = data[180:359,0:nl-1]
; out_arr[fi,180:359,1:nl] = data[0:179,0:nl-1]
; end
; 'MRI-ESM2-0': begin
; out_arr[fi, 0:179,0:nl-1] = data[180:359,0:nl-1]
; out_arr[fi,180:359,0:nl-1] = data[0:179,0:nl-1]
; end
'NorESM2-LM': begin
out_arr[fi,0:178,0:nl-2] = data[179:357,0:nl-2] ;[0,90],[358,181]
out_arr[fi,179,0:nl-2] = data[0,0:nl-2]
out_arr[fi,180:358,0:nl-2] = data[0:178,0:nl-2]
out_arr[fi,359,0:nl-2] = data[179,0:nl-2]
end
'BCC-CSM2-MR': begin
out_arr[fi,0:179,1:nl] = data[179:358,0:nl-1] ;[0,89],[359,179]
out_arr[fi,180:359,1:nl] = data[0:179,0:nl-1]
end
'CMCC-ESM2': begin
out_arr[fi,0:179,0:nl-2] = data[179:358,0:nl-2] ;[0,90],[359,181]
out_arr[fi,180:359,0:nl-2] = data[0:179,0:nl-2]
end
'ACCESS-CM2': begin
out_arr[fi,0:179,1:nl] = data[179:358,0:nl-1] ;[1,89],[359,179]
out_arr[fi,180:359,1:nl] = data[0:179,0:nl-1]
end
'CanESM5': begin
out_arr[fi,0:179,2:nl+1] = data[178:357,0:nl-1] ;[0,88],[358,176]
out_arr[fi,180:359,2:nl+1] = data[0:179,0:nl-1]
end
'CESM2-WACCM': begin
out_arr[fi,0:179,0:nl-2] = data[179:358,0:nl-2] ;[0,90],[359,181]
out_arr[fi,180:359,0:nl-2] = data[0:179,0:nl-2]
end
'CMCC-CM2-SR5': begin
out_arr[fi,0:179,0:nl-2] = data[179:358,0:nl-2] ;[0,90],[359,181]
out_arr[fi,180:359,0:nl-2] = data[0:179,0:nl-2]
end
'FGOALS-g3': begin
out_arr[fi,0:179,0:nl-2] = data[179:358,0:nl-2] ;[0,90],[359,181]
out_arr[fi,180:359,0:nl-2] = data[0:179,0:nl-2]
end
'MIROC6': begin
out_arr[fi,0:179,1:nl] = data[179:358,0:nl-1] ;[0,89],[359,178]
out_arr[fi,180:359,1:nl] = data[0:179,0:nl-1]
end
'MPI-ESM1-2-HR': begin
out_arr[fi,0:179,1:nl] = data[180:359,0:nl-1] ;[0,89],[360,179]
out_arr[fi,180:359,1:nl] = data[0:179,0:nl-1]
end
'MPI-ESM1-2-LR': begin
out_arr[fi,0:179,1:nl] = data[179:358,0:nl-1] ;[0,89],[359,178]
out_arr[fi,180:359,1:nl] = data[0:179,0:nl-1]
end
'MRI-ESM2-0': begin
out_arr[fi,0:179,1:nl] = data[179:358,0:nl-1] ;[0,89],[359,179]
out_arr[fi,180:359,1:nl] = data[0:179,0:nl-1]
end
endcase
year = start_year + fi/12;
month= fi mod 12 +1.0
date = year*100.+month
fname = perfname+string(date,format='(g0)')+'.tif'
out_path = writepath1+fname
WRITE_TIFF,out_path,round(out_arr[fi,*,*]*10000.0)/10000.0,GEOTIFF=GEOTIFF,/FLOAT
endfor
;依次输出年份
writepath2 = writepath+'annual\'
if file_test(writepath2, /directory) eq 0 then file_mkdir,writepath2;输出路径没有即新建立
for fi = 0,nb_data/12-1 do begin
;print,fi
year = start_year + fi;
out_tiff = fltarr(360,180)
ns=N_elements(out_tiff[*,0]);获取列数
nl=N_elements(out_tiff[0,*]);获取行数
for i=0,ns-1 do begin
for j=0,nl-1 do begin
out_tiff[i,j] = mean(out_arr[fi:fi+11,i,j])
if out_tiff[i,j] lt -100 then out_tiff[i,j] = -999
endfor
endfor
date = year
fname = perfname+string(date,format='(g0)')+'.tif'
out_path = writepath2+fname
WRITE_TIFF,out_path,round(out_tiff*10000.0)/10000.0,GEOTIFF=GEOTIFF,/FLOAT
endfor
;;;;;;;释放内存
;fids = envi_get_file_ids();获取当前内存中的所有文件的fid
;size = N_elements(fids);获取数组的大小
;length = size
;;循环释放内存中的文件
;for i = 0L, length-1 do begin
; envi_file_mng,id = fids[i],/remove
;endfor
temp=0;
end
