====================================================================== ; ESMF_regrid_35.ncl ; ; Concepts illustrated: ; - Interpolating from high-to-low resolution using 3 ESMF_regrid methods ; - Determine the number of _FillValue (nmsg); use 'mask' if nmsg>0 ; - Panel plot ; - Create netCDF ;====================================================================== ; ; These files are loaded by default in NCL V6.2.0 and newer ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ; ; This file still has to be loaded manually load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl" begin PLOT = True NC = True ;---Read data off source file hfcVarName = "E_HFC_1234" ; variable name hfcDirName = "./" ; source file directory hfcFileName = "HFC.nc" ; source grid hfcPathName = hfcDirName+hfcFileName hfile = addfile(hfcPathName,"r") x = hfile->$hfcVarName$ printVarSummary(x) printMinMax(x,0) nmsg = num(ismissing(x)) ; Check to see id any _FillValue print("nmsg="+nmsg) ;---Look at statistics; help with plotting range opt = True opt@PrintStat = True statb = stat_dispersion(x, opt ) ;---Create destination 1x1 grid ; latitude: [33.04999923706055..72.94999694824219] ; longitude: [-12.45001220703125..44.95000076293945] nlat = 41 lat = fspan(33, 73, nlat) lat@units = "degrees_north" lat!0 = "lat" lat&lat = lat mlon = 58 lon = fspan(-12, 45, mlon) lon@units = "degrees_north" lon!0 = "lon" lon&lon = lon ;---Assign zoom region for graphics minlon = min(x&longitude) maxlon = max(x&longitude) minlat = min(x&latitude) maxlat = max(x&latitude) print("min/max lat = " + minlat + "/" + maxlat) print("min/max lon = " + minlon + "/" + maxlon) ;---Options for regridding Opt = True Opt@ForceOverwrite = True Opt@SrcTitle = hfcFileName ; source grid if (nmsg.gt.0) then Opt@SrcMask2D = where(ismissing(x),0,1) end if Opt@SrcRegional = True Opt@DstTitle = "HFC 1x1 Rectilinear Grid" ; destination grid Opt@DstGridLat = lat Opt@DstGridLon = lon Opt@DstRegional = True ;;Opt@PrintTimings = True ;---Regrid Opt@InterpMethod = "bilinear" ; bilinear interpolation Opt@WgtFileName = "HFC_1x1_bilinear.nc" x_regrid_b = ESMF_regrid(x,Opt) Opt@InterpMethod = "patch" ; patch interpolation Opt@WgtFileName = "HFC_1x1_patch.nc" x_regrid_p = ESMF_regrid(x,Opt) Opt@InterpMethod = "conserve" ; interpolation Opt@WgtFileName = "HFC_1x1_conserve.nc" x_regrid_c = ESMF_regrid(x,Opt) ;---------------------------------------------------------------------- ; Plotting section ;---------------------------------------------------------------------- if (PLOT) then wks = gsn_open_wks("png","ESMF_regrid") ; send graphics to PNG file res = True ; Plot mods desired. res@gsnDraw = False res@gsnFrame = False res@mpDataBaseVersion = "MediumRes" res@mpDataSetName = "Earth..4" res@mpFillOn = False res@mpOutlineOn = True res@mpOutlineBoundarySets = "AllBoundaries" res@mpMinLatF = min(lat) res@mpMaxLatF = max(lat) res@mpMinLonF = min(lon) res@mpMaxLonF = max(lon) res@cnLinesOn = False res@cnFillMode = "RasterFill" res@cnLineLabelsOn = False res@cnFillOn = True res@cnFillPalette = "WhiteBlueGreenYellowRed" res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = (/0.001, 0.050, 0.060, 0.070, 0.080, 0.090, 0.100, 0.125, 0.150\ ,0.200, 0.250, 0.50, 1, 5, 10, 15, 20, 25, 50 /) res@pmTickMarkDisplayMode = "Always" res@gsnAddCyclic = False ; don't add cyclic longitude point res@tiMainFontHeightF = 0.02 res@lbLabelBarOn = False res@gsnLeftString = "Original data (" + \ str_join(tostring(dimsizes(x))," x ") + ")" plot = gsn_csm_contour_map(wks,x,res) ;---bilinear res@gsnLeftString = "Regridded data (bilinear) (" + \ str_join(tostring(dimsizes(x_regrid_b))," x ") + ")" plot_b = gsn_csm_contour_map(wks,x_regrid_b,res) ;---patch res@gsnLeftString = "Regridded data (patch) (" + \ str_join(tostring(dimsizes(x_regrid_p))," x ") + ")" plot_p = gsn_csm_contour_map(wks,x_regrid_p,res) ;---conserve res@gsnLeftString = "Regridded data (conserve) (" + \ str_join(tostring(dimsizes(x_regrid_c))," x ") + ")" plot_c = gsn_csm_contour_map(wks,x_regrid_c,res) ;---Panel all four plots pres = True pres@gsnMaximize = True pres@gsnPanelLabelBar = True pres@pmLabelBarWidthF = 0.8 pres@lbLabelFontHeightF = 0.01 gsn_panel(wks,(/plot,plot_b,plot_p,plot_c/),(/2,2/),pres) end if ; PLOT if (NC) then ; Create netCDF diro = "./" filo = "HFC_ESMF_REGRID.nc" system("/bin/rm -f "+diro+filo) ; remove any pre-existing file ncdf = addfile(diro+filo ,"c") ; open output netCDF file ncdf@title = "HFC: ESMF Regrid" ncdf@source = hfcFileName ncdf@creation_date = systemfunc("date") name = hfcVarName+"_bilinear" ncdf->$name$ = x_regrid_b name = hfcVarName+"_patch" ncdf->$name$ = x_regrid_p name = hfcVarName+"_conserve" ncdf->$name$ = x_regrid_c end if end