;---Specify variable [%> ncl_filedump .....] vName = "APCP_P8_2L1_GLC0_acc6h" ;---Specify interpolation method InterpMethod = "bilinear" ; "patch", "conserve", "bilinear" ; "neareststod" ,"nearestdtos" ;---Write regridded variable to netCDF netCDF= True if (netCDF) then ; the following are ignored if netCDF=False rgrdDir = "./" rgrdFile = "QPF_regrid."+InterpMethod+".nc" rgrdPath = rgrdDir + rgrdFile end if ;---Specify graphics file stuff pltType = "png" pltDir = "./" pltName = "QPF."+vName+"."+InterpMethod pltPath = pltDir+pltName ;-- File that contains data and lat/lon source grid srcDir = "./" srcFileName = "sourceGrid_stovern.grb" srcPathName = srcDir+srcFileName ;---File that contains lat/lon destination grid dstDir = "./" dstFileName = "destinationGrid_stovern.grb" dstPathName = dstDir+dstFileName ;---Open both files sfile = addfile(srcPathName,"r") dfile = addfile(dstPathName,"r") ;---Get variable off source file var = sfile->$vName$ ; also: sfile->APCP_P8_2L1_GLC0_acc6h var_rank = dimsizes(dimsizes(var)) ;---Read source lat/lon values src_lat2d = sfile->gridlat_0 src_lon2d = sfile->gridlon_0 ;---Read destination lat/lon values dst_lat2d = dfile->gridlat_0 dst_lon2d = dfile->gridlon_0 ;---Dimension [size] information src_dim = dimsizes(src_lat2d) src_nlat = src_dim(0) src_nlon = src_dim(1) dst_dim = dimsizes(dst_lat2d) dst_nlat = dst_dim(0) dst_nlon = dst_dim(1) ;---Map limits: ; Used for plotting later dst_minlat = min(dst_lat2d) dst_maxlat = max(dst_lat2d) dst_minlon = min(dst_lon2d) dst_maxlon = max(dst_lon2d) ;---Debug prints print("-----src_lat2d----") print("min/max src_lat2d = " + min(src_lat2d) + "/" + max(src_lat2d)) print("min/max src_lon2d = " + min(src_lon2d) + "/" + max(src_lon2d)) print("min/max dst_lat2d = " + dst_minlat + "/" + dst_maxlat) print("min/max dst_lon2d = " + dst_minlon + "/" + dst_maxlon) ;---Set up options for regridding Opt = True Opt@WgtFileName = "QPF_weight_"+InterpMethod+".nc" Opt@SrcGridLat = src_lat2d ; source grid Opt@SrcGridLon = src_lon2d Opt@DstGridLat = dst_lat2d ; destination grid Opt@DstGridLon = dst_lon2d Opt@SrcRegional = True ; Necessary if grids Opt@DstRegional = True ; are regional Opt@InterpMethod = InterpMethod Opt@RemoveSrcFile = True Opt@RemoveDstFile = True Opt@ForceOverwrite = True ; Optional, but recommended. Opt@PrintTimings = True ; Optional. Opt@Debug = True ; Optional ;---ESMF_regrid is an externally developed package that is invoked by NCL wcBeginRegrid = systemfunc("date") var_regrid = ESMF_regrid(var,Opt) ; Do the regridding wallClockElapseTime(wcBeginRegrid, "===>ESMF_regrid time: ", 0) print("") ;---More debug prints print("---------var----------") printMinMax(var,0) print("-------var_regrid---------") printVarSummary(var_regrid) printMinMax(var_regrid,0) ;---------------------------------------------------------------------- ; Plot the original and regridded data in a panel plot. ;---------------------------------------------------------------------- var@lat2d = src_lat2d ; source coordinates var@lon2d = src_lon2d ;---Used in plot title below src_dims = tostring(src_dim) dst_dims = tostring(dst_dim) ;---For 'fun' look at the distribution of the data. Can be useful for graphics opt_sd = True opt_sd@PrintStat = True var_stat = stat_dispersion(var, opt_sd ) begPlot = get_cpu_time() wks = gsn_open_wks(pltType,pltPath) ; send graphics to PNG file res = True ; Plot mods desired. ;---General resources res@gsnDraw = False ; We will panel later. res@gsnFrame = False res@gsnMaximize = True ; Maximize plot res@cnFillMode = "AreaFill" ; Default ;res@cnFillMode = "RasterFill" ; Raster Mode ;res@cnFillMode = "CellFill" ; Cell Mode here faster than Raster ;---Color Palette ; http://www.ncl.ucar.edu/Document/Graphics/color_table_gallery.shtml#Starts-with-white res@cnFillPalette = "prcp_1" res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = (/0.1,0.5,1.0,1.5,2.0,3,4,5,7.5,10,20,30,40,50/) ;---Contour and labelbar resources res@cnFillOn = True ; Color plot desired. res@cnLinesOn = False ; Turn off contour lines. res@cnLineLabelsOn = False ; Turn off contour labels. res@lbLabelBarOn = False ; Labelbar will be in panel. res@gsnAddCyclic = False ; Data is regional. ;---Map resources res@mpFillOn = False ;res@mpDataBaseVersion = "MediumRes" ;res@mpOutlineBoundarySets = "AllBoundaries" ;res@mpUSStateLineColor = "gray" ;res@mpGeophysicalLineColor = "gray" ;res@mpNationalLineColor = "gray" ;---Zoom in on map mp_offset = 0.03 res@mpMinLatF = dst_minlat-mp_offset ; Leave a little bit res@mpMinLonF = dst_minlon-mp_offset ; of a margin. res@mpMaxLatF = dst_maxlat+mp_offset res@mpMaxLonF = dst_maxlon+mp_offset ;---Titles and tickmarks res@pmTickMarkDisplayMode = True ; Better tickmarks res@tmYROn = False res@tmXTOn = False res@tmXBLabelFontHeightF = 0.01 res@tmYLLabelFontHeightF = 0.01 res@gsnLeftString = "" res@gsnRightString = "" res@tiMainFontHeightF = 0.020 ;res@tiMainOffsetYF = -0.02 ; move closer to plot ;res@gsnStringFontHeightF = 0.015 ;---Resources for paneling pres = True pres@gsnMaximize = True pres@gsnPanelLabelBar = True pres@pmLabelBarWidthF = 0.8 ; SAME_COLORBAR = False if(SAME_COLORBAR) then mnmxint = nice_mnmxintvl( min(var_regrid), max(var_regrid), 18, False) res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = mnmxint(0) res@cnMaxLevelValF = mnmxint(1) res@cnLevelSpacingF = mnmxint(2) end if res@tiMainString = "Regridded data (" + \ str_join(dst_dims," x ") + ")" plot_regrid = gsn_csm_contour_map(wks,var_regrid,res) ;---Make sure same contours levels are used for both plots. if(.not.SAME_COLORBAR) then getvalues plot_regrid@contour "cnLevels" : levels end getvalues res@cnLevels = levels res@cnLevelSelectionMode = "ExplicitLevels" end if res@tiMainString = "Original data (" + \ str_join(src_dims," x ") + ")" plot_orig = gsn_csm_contour_map(wks,var,res) pres@gsnPanelMainString = var@long_name +"[ "+var@units+"]: "+ var@initial_time+ ": ftime="+var@forecast_time gsn_panel(wks,(/plot_orig,plot_regrid/),(/2,1/),pres) print("Plot time "+res@cnFillMode+": " + (get_cpu_time() - begPlot)) print("") ;---Write to netCDF if (netCDF) then delete([/var_regrid@lat2d, var_regrid@lon2d /]) system("rm -f " + rgrdPath) ; remove any preexisting file rgrd_nc = addfile(rgrdPath,"c") ;---Create variable to hold global file attributes global = True global@regrid = "NCL: ESMF_regrid(NCL version '" + \ get_ncl_version() + "')" global@regrid_method = InterpMethod global@creation_date = systemfunc("date") fileattdef( rgrd_nc, global ) ; copy global file attributes ;;;;filedimdef(rgrd_nc,"TIME",-1,True) ; force an unlimited dimension rgrd_nc->$vName$ = var_regrid rgrd_nc->lat2d = dst_lat2d rgrd_nc->lon2d = dst_lon2d end if