pro CMIP6_batch_pr

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']

;me

;model_name = 'NorESM2-LM'?

;var_name = ['pr','lon','lat']

;model_name = 'BCC-CSM2-MR'

;var_name = ['pr','lon','lat']

;

;model_name = 'CMCC-ESM2'

;var_name = ['pr','lon','lat']

;

;model_name = 'ACCESS-CM2'

;var_name = ['pr','lon','lat']

;model_name = 'CanESM5'

;var_name = ['pr','lon','lat']

;model_name = 'CESM2-WACCM'

;var_name = ['pr','lon','lat']

;model_name = 'CMCC-CM2-SR5'

;var_name = ['pr','lon','lat']?

;model_name = 'FGOALS-g3'

;var_name = ['pr','lon','lat']?

;model_name = 'MIROC6'

;var_name = ['pr','lon','lat']?

;model_name = 'MPI-ESM1-2-HR'

;var_name = ['pr','lon','lat']

;model_name = 'MPI-ESM1-2-LR'

;var_name = ['pr','lon','lat']?

model_name = 'MRI-ESM2-0'

var_name = ['pr','lon','lat']

file_arr = file_search(netcdf_path,'pr*'+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]+'_'+'pr_'

? 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_pr\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參數(shù)：顯示處理進(jìn)度條

? if file_test(glt_path, /directory) eq 0 then file_mkdir,glt_path;輸出路徑?jīng)]有即新建立

? 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;? 關(guān)閉

? temp= 0 ;

??

? nt=N_elements(tos[0,0,*]);獲取時(shí)間數(shù)

? ns=N_elements(tos[*,0,0]);獲取列數(shù)

? nl=N_elements(tos[0,*,0]);獲取行數(shù)

??

? 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

? ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;存儲(chǔ)經(jīng)緯度數(shù)據(jù)并構(gòu)建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

? ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;存儲(chǔ)tos數(shù)據(jù)并做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;

;;;;;;釋放內(nèi)存

fids = envi_get_file_ids();獲取當(dāng)前內(nèi)存中的所有文件的fid

size = size(fids);獲取數(shù)組的大小

length = size[1]

;循環(huán)釋放內(nèi)存中的文件

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

;;;;;所有的數(shù)據(jù)統(tǒng)一為-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;輸出路徑?jīng)]有即新建立??

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]);獲取列數(shù)

? nl=N_elements( data[0,*]);獲取行數(shù)

? 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 0 then data[i,j] = -999

? endfor

? endfor

??

? 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;輸出路徑?jīng)]有即新建立

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]);獲取列數(shù)

? nl=N_elements(out_tiff[0,*]);獲取行數(shù)

? 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 0 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??

??

;;;;;;;釋放內(nèi)存

;fids = envi_get_file_ids();獲取當(dāng)前內(nèi)存中的所有文件的fid

;size = N_elements(fids);獲取數(shù)組的大小

;length = size

;;循環(huán)釋放內(nèi)存中的文件

;for i = 0L, length-1 do begin

;? envi_file_mng,id = fids[i],/remove

;endfor

??

??

? temp=0;

??

end