网格气象场插值-NCL版

load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
 ;-- load ESMF library
begin
f = addfile("rectilinear_grid_2D.nc","r") ;-- open file
var = f->tsurf(0,:,:) ;-- read variable
;-- define the T255 grid (lat x lon - 256x512)
nlat = 256 ;-- number of latitudes
nlon = 512 ;-- number of longitudes

grint = 0.54 ;-- grid spacing
dst_lat = fspan((-90.0+grint),(90.-grint),nlat)*1d ;-- type double
dst_lat!0 = "lat" ;-- dimension name
dst_lat@units = "degrees_north" ;-- dimension units
dst_lon = fspan(0.0,(360.-grint),nlon)*1d ;-- type double
dst_lon!0 = "lon" ;-- dimension name
dst_lon@units = "degrees_east" ;-- dimension units

;-- set ESMF resources
Opt = True
Opt@InterpMethod = "bilinear" ;-- interpolation method
Opt@SrcFileName = "ECHAM5_SCRIP_bilinear.nc" ;-- new source file name
Opt@WgtFileName = "ECHAM5toWorldCurvilinear_bilinear.nc";-- weights file
Opt@ForceOverwrite = True ;-- force overwrite
Opt@DstFileName = "WorldRectilinear_SCRIP_bilinear.nc";-- destination file
Opt@DstGridType = "rectilinear" ;-- Destination grid
Opt@DstGridLon = dst_lon
Opt@DstGridLat = dst_lat
var_regrid = ESMF_regrid(var,Opt) ;-- regrid the variable
var_regrid!0 = "y" ;-- name dimension 0 (default: ncl0)
var_regrid!1 = "x" ;-- name dimension 1 (default: ncl1)

ofile = "regrid_T63_to_T255.nc"
system("rm -f "+ofile) ;-- delete netCDF file if it exist
cdf_file = addfile(ofile,"c") ;-- create a new netCDF file
delete_VarAtts(var_regrid,(/"lat2d", "lon2d"/)) ;-- delete the attributes
 ;-- lat2d and lon2d
cdf_file->lat2d = lat2d ;-- write lat2d to file
cdf_file->lon2d = lon2d ;-- write lon2d to file
cdf_file->S = var_regrid ;-- write variable to file
end

load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
begin
fili = "rectilinear_grid_2D.nc" ;-- data file name
f = addfile(fili,"r") ;-- open data file
var = f->tsurf(0,:,:) ;-- read variable

;-- read the destination grid lat/lon arrays
dstfili = "curvilinear_ocean.nc" ;-- destination grid file
d = addfile(diri+dstfili,"r") ;-- open grid file
dst_lat = d->lat ;-- 2D latitudes
dst_lon = d->lon ;-- 2D longitudes
dims = dimsizes(dst_lat) ;-- size of lat/lon
nlat = dims(0) ;-- number of latitudes
nlon = dims(1) ;-- number of longitudes

;-- set ESMF resources
Opt = True
Opt@InterpMethod = "bilinear" ;-- interpolation method
Opt@SrcFileName = "T63_SCRIP_bilinear.nc" ;-- source file name
Opt@DstFileName = "WorldCurvilin_SCRIP_bilin.nc" ;-- destination file
Opt@WgtFileName = "T63toWorldCurvilinear_bilinear.nc"
 ;-- name of weights file, which will will be generated
Opt@ForceOverwrite = True ;-- force overwrite
Opt@DstMask2D = where(ismissing(dvar),0,1) ;-- create mask from
 ;-- destination
Opt@DstGridType = "curvilinear" ;-- destination grid
Opt@DstTitle = "World Grid Curvilinear Resolution bilinear"
 ;-- destination title
Opt@DstGridLon = dst_lon ;-- set destination lon
Opt@DstGridLat = dst_lat ;-- set destination lat
;-- call ESMF_regrid
var_regrid = ESMF_regrid(var,Opt) ;-- regrid variable
var_regrid!0 = "y" ;-- named coordinate dimension 0
var_regrid!1 = "x" ;-- named coordinate dimension 1
 delete_VarAtts(var_regrid,(/"lat2d", "lon2d"/)) ;-- delete attributes
;-- assign a output netcdf file for the new regridded data name of output file
ofile = "regridd_rectilin_to_curvilinear_bilin_wgts_destgrid_ESMF.nc"
system("rm -rf "+ofile) ;-- delete ofile if it exists
fout = addfile(ofile, "c") ;-- open new file

;-- start to define output file settings
setfileoption(fout,"DefineMode",True) ;-- explicitly declare file def mode
;-- create global attributes of output file
fAtt = True ;-- assign file attributes
fAtt@Conventions = "CF-1.5"
fAtt@title = "Regrid T63 to curvilinear grid (MPIOM)"
fAtt@source_file = fili
fAtt@creation_date = systemfunc ("date")
fAtt@history = \
 "ncl NCL_Advanced_regrid_rectilin_to_curvilin_bilin_wgts_ESMF.ncl"
fileattdef(fout,fAtt) ;-- copy file attributes

;-- predefine the coordinate variables and their dimensionality
dimNames = (/"y", "x"/) ;-- curvilinear grid: dimensions not lat/lon
dimSizes = (/nlat, nlon/) ;-- dimension size of destination y/x
dimUnlim = (/False, False/)
filedimdef(fout,dimNames,dimSizes,dimUnlim) ;-- define dimensions

;-- predefine the the dimensionality of the variables to be written out
filevardef(fout,"lat",typeof(dst_lat),getvardims(dst_lat)) ;-- variable lat
filevardef(fout,"lon",typeof(dst_lon),getvardims(dst_lon)) ;-- variable lon
filevardef(fout,"var",typeof(var_regrid),getvardims(var_regrid))
;-- copy attributes associated with var to output file
;-- copy variable attributes
filevarattdef(fout,"lat", dst_lat) ;-- copy attributes from destination lat
 filevarattdef(fout,"lon", dst_lon) ;-- copy attributes from destination lon
 filevarattdef(fout,"var", var_regrid) ;-- copy var_regrid attributes
;-- explicitly exit file definition mode (not required)
setfileoption(fout,"DefineMode",False)
;-- output only the data values since the dimensionality and such have been
;-- predefined; the "(/ /)" syntax tells NCL to only output the data values
fout->var = (/var_regrid/) ;-- write variable to new netCDF file
fout->lat = (/dst_lat/) ;-- write lat to new netCDF file
fout->lon = (/dst_lon/) ;-- write lon to new netCDF file

end

load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
begin
fili = "curvilinear_ocean.nc"
 sfile = addfile(fili,"r") ;-- open data file
thetao = sfile->thetao(0,0,:,:) ;-- read variable
thetao@lat2d = sfile->lat ;-- tells ESMF what kind of source
 ;-- grid
 thetao@lon2d = sfile->lon ;-- tells ESMF what kind of source
 ;-- grid
;-- set resources
Opt = True
Opt@InterpMethod = "bilinear" ;-- interpolation method
Opt@SrcFileName = "Curvilin_SCRIP_bilinear.nc" ;-- source file name
Opt@DstFileName = "World1deg_SCRIP_bilinear.nc" ;-- destination file
Opt@WgtFileName = "CurvilintoWORLD_1x1_bilinear.nc"
 ;-- name of weights file, which will be generated
Opt@ForceOverwrite = True ;-- force overwrite

Opt@SrcMask2D = where(.not. ismissing(thetao),1,0) ;-- if data contains
 ;-- missing values
Opt@DstGridType = "1deg" ;-- Destination grid
Opt@DstTitle = "World Grid 1x1-degree Resolution bilinear"
 ;-- destination title
Opt@DstLLCorner = (/-89.5d, 0.0d /) ;-- destination lower
 ;-- left corner
Opt@DstURCorner = (/ 89.5d, 359.0d /) ;-- destination upper
 ;-- right corner
;-- call ESMF_regrid
thetao_regrid = ESMF_regrid(thetao,Opt)
nlon = dimsizes(thetao_regrid&lon)
nlat = dimsizes(thetao_regrid&lat)
;-- assign a output netcdf file for the new regridded data;
;-- (npoints = 180x360)
ofile = "regridded_curvilinear_bilinear_1x1deg_thetao_ESMF.nc"
system("rm -rf "+ofile) ;-- delete file if exists
fout = addfile(ofile,"c") ;-- open new file
;-- start to define output file settings
setfileoption(fout,"DefineMode",True)
 ;-- explicitly declare file definition mode

 ;-- create global attributes of the file
fAtt = True ;-- assign file attributes
fAtt@Conventions = "CF-1.5“
fAtt@project_id = "DKRZ NCL Advanced Workshop – Regridding"
fAtt@source_file = fili
fAtt@creation_date = systemfunc ("date")
fAtt@history = \
 "ncl NCL_Advanced_regrid_curvilin_to_rectilin_bilin_weights_ESMF.ncl"
fileattdef(fout,fAtt) ;-- copy file attributes
;-- predefine the coordinate variables and their dimensionality
dimNames = (/"lat", "lon"/)
dimSizes = (/nlat, nlon/)
dimUnlim = (/False, False/)
filedimdef(fout,dimNames,dimSizes,dimUnlim) ;-- define dimensions
;-- predefine the the dimensionality of the variables to be written out
filevardef(fout,"lat",typeof(thetao_regrid&lat),getvardims(thetao_regrid&lat))
filevardef(fout,"lon",typeof(thetao_regrid&lon),getvardims(thetao_regrid&lon))
filevardef(fout,"thetao",typeof(thetao_regrid),getvardims(thetao_regrid))
;-- copy attributes associated with each variable to the file
filevarattdef(fout,"lat",thetao_regrid&lat) ;-- copy lat attributes
filevarattdef(fout,"lon",thetao_regrid&lon) ;-- copy lon attributes
filevarattdef(fout,"thetao",thetao_regrid)
 ;-- copy thetao_regrid attributes

 ;-- explicitly exit file definition mode (not required)
setfileoption(fout,"DefineMode",False)
;-- output only the data values since the dimensionality and such have been
;-- predefined; the "(/ /)" syntax tells NCL to only output the data values
fout->thetao = (/thetao_regrid/) ;-- write variable to new netCDF file
fout->lat = (/thetao_regrid&lat/) ;-- write lat to new netCDF file
fout->lon = (/thetao_regrid&lon/) ;-- write lon to new netCDF file

end

load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
begin
rad2deg = 45./atan(1.) ;-- radians to degrees
f = addfile("triangular_grid_ICON.nc","r")
var = f->S(time|0,depth|0,ncells|:) ;-- read variable
x = f->clon * rad2deg ;-- cell center, lon
 y = f->clat * rad2deg ;-- cell center, lat
x!0 = "lon" ;-- set named dimension lon
y!0 = "lat" ;-- set named dimension lat
x@units = "degrees_east" ;-- set lon units
y@units = "degrees_north" ;-- set lat units
;-- set ESMF resources
Opt = True
Opt@InterpMethod = "bilinear" ;-- interpolation method
Opt@ForceOverwrite = True ;-- force overwrite
Opt@SrcFileName = "Unstruct_SCRIP_bilinear.nc";-- source file name 
Opt@SrcInputFileName = diri+fili ;-- optional, but good idea
Opt@SrcRegional = False
Opt@SrcGridLat = y
Opt@SrcGridLon = x
Opt@WgtFileName = "UnstructtoWORLD_1x1_bilinear.nc"
 ;-- name of weights file
 Opt@DstFileName = "World1deg_SCRIP_bilinear.nc" ;-- destination file name
Opt@DstGridType = "rectilinear" ;-- destination grid
Opt@DstTitle = "World Grid 1x1-degree Resolution bilinear"
 ;-- dest. title
Opt@DstRegional = False
Opt@DstGridLon = fspan(-180.,180.,360)
Opt@DstGridLat = fspan(-90.,90.,180) ;-- call ESMF_regrid
var_regrid = ESMF_regrid(var,Opt) ;-- do the regridding
nlon = dimsizes(var_regrid&lon) ;-- dim size new lon
nlat = dimsizes(var_regrid&lat) ;-- dim size new lat
;-- assign a output netcdf file for the new regridded data (npoints = 180x360)
system("rm -rf regridded_unstructured_to_rectilinear_bilinear_ESMF.nc")
fout = addfile("regrid_unstructured_to_rectilin_bilin_ESMF.nc", "c")
;-- start to define output file settings
setfileoption(fout,"DefineMode",True) ;-- explicitly declare file def. mode

;-- create global attributes of the file
fAtt = True ;-- assign file attributes
fAtt@Conventions = "CF-1.5"
fAtt@source_file = fili
fAtt@creation_date = systemfunc ("date")
fAtt@history = \
 "ncl NCL_Advanced_regrid_unstruct_to_rectilin_bilinear_wgts_ESMF.ncl"
fileattdef(fout,fAtt) ;-- copy file attributes
;-- predefine the coordinate variables and their dimensionality
dimNames = (/"lat", "lon"/)
dimSizes = (/nlat, nlon/)
dimUnlim = (/False, False/)
filedimdef(fout,dimNames,dimSizes,dimUnlim)
;-- predefine the the dimensionality of the variables to be written out
filevardef(fout, "lat", typeof(var_regrid&lat), getvardims(var_regrid&lat))
filevardef(fout, "lon", typeof(var_regrid&lon),getvardims(var_regrid&lon))
 filevardef(fout, "S", typeof(var_regrid), getvardims(var_regrid))

 ;-- copy attributes associated with each variable to the file
filevarattdef(fout,"lat", var_regrid&lat) ;-- copy lat attributes
filevarattdef(fout,"lon", var_regrid&lon) ;-- copy lon attributes
filevarattdef(fout,"S", var_regrid) ;-- copy var_regrid attributes
;-- explicitly exit file definition mode (not required)
setfileoption(fout,"DefineMode",False)
;-- output only the data values since the dimensionality and such have
;-- been predefined; the "(/ /)" syntax tells NCL to only output the
;-- data values to the predefined locations on the file.
fout->S = (/var_regrid/) ;-- write variable to new netCDF file
fout->lat = (/var_regrid&lat/) ;-- write lat to new netCDF file
fout->lon = (/var_regrid&lon/) ;-- write lon to new netCDF file
end