wrfout输出数据使用ncl进行插值到ERA_Interim数据上，数据与地图的匹配问题

王精志

begin
path1 = "/Volumes/UNTITLED/1979/"
p1    = systemfunc("cd " + path1 +"; ls wrfout_d01_1979-0*")
;print(p1)
d = dimsizes(p1)
f1    = addfiles(path1+p1, "r")
ListSetType(f1, "cat")
printVarSummary(f1)

f = addfile ("/Volumes/UNTITLED/1979/1979.nc", "r")
f2      = addfile ("/Volumes/UNTITLED/1979/sfc.to.nc", "r")
LAT= wrf_user_getvar(f1,"XLAT",0)   ; latitude/longitude
LON = wrf_user_getvar(f1,"XLONG",0)  ; required for plotting

;lon = f->longitude({min(LON):max(LON)})
;lat = f->latitude({min(LAT):max(LAT)})

longitude = lonFlip(f->longitude)
latitude = f->latitude

lon = longitude({min(LON):max(LON)})
lat = latitude({min(LAT):max(LAT)})

;土壤温度
TSLB = wrf_user_getvar(f1,"TSLB",-1)
tslb_xgrd  = rcm2rgrid_Wrap(LAT,LON,TSLB,lat,lon,0)

tslb_xgrd@lat2d = LAT
tslb_xgrd@lon2d = LON
printVarSummary(TSLB)
tslb1 = tslb_xgrd(:,0,:,:)
tslb2 = tslb_xgrd(:,1,:,:)
tslb3 = tslb_xgrd(:,2,:,:)
tslb4 = tslb_xgrd(:,3,:,:)

printVarSummary(tslb1)

;============================================================================================================================
wks1 = gsn_open_wks("png","/Volumes/UNTITLED/1979/spin_up_1979(TSLB)_map")

;tslb1_xgrd  = rcm2rgrid_Wrap(LAT,LON,tslb1,lat,lon,0)
;tslb2_xgrd  = rcm2rgrid_Wrap(LAT,LON,tslb2,lat,lon,0)
;tslb3_xgrd  = rcm2rgrid_Wrap(LAT,LON,tslb3,lat,lon,0)
;tslb4_xgrd  = rcm2rgrid_Wrap(LAT,LON,tslb4,lat,lon,0)

;delete(tslb1_xgrd&latitude)
;delete(tslb1_xgrd&longitude)

tslb1@lat2d = LAT
tslb1@lon2d = LON

daily_avg_tslb1_t = new((/4,55,98/),"float")
daily_avg_tslb1_t(0,:,:) = dim_avg_n_Wrap(tslb1,(/0/))
daily_avg_tslb1_t(1,:,:) = dim_avg_n_Wrap(tslb2,(/0/))
daily_avg_tslb1_t(2,:,:) = dim_avg_n_Wrap(tslb3,(/0/))
daily_avg_tslb1_t(3,:,:) = dim_avg_n_Wrap(tslb4,(/0/))

printVarSummary(daily_avg_tslb1_t)

a = daily_avg_tslb1_t(0,:,:)
delete(a@lat2d)
delete(a@lon2d)
printVarSummary(LAT)
stl1 = short2flt(f->stl1(:,{min(LAT):max(LAT)},{min(LON):max(LON)}))
printVarSummary(stl1)
stl2 = short2flt(f->stl2(:,{min(LAT):max(LAT)},{min(LON):max(LON)}))
stl3 = short2flt(f->stl3(:,{min(LAT):max(LAT)},{min(LON):max(LON)}))
stl4 = short2flt(f->stl4(:,{min(LAT):max(LAT)},{min(LON):max(LON)}))

daily_avg_tslb1_t_ERA = new((/4,55,98/),"float")
daily_avg_tslb1_t_ERA(0,:,:) = dim_avg_n_Wrap(stl1,(/0/))
daily_avg_tslb1_t_ERA(1,:,:) = dim_avg_n_Wrap(stl2,(/0/))
daily_avg_tslb1_t_ERA(2,:,:) = dim_avg_n_Wrap(stl3,(/0/))
daily_avg_tslb1_t_ERA(3,:,:) = dim_avg_n_Wrap(stl4,(/0/))


res               = True
res@gsnMaximize   = True   ; maximize plot in frame
res@gsnFrame          =       False
res@gsnDraw           =       False
res@cnFillOn            =True        ; turn on color

res@cnInfoLabelOn        = False
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels = (/264,268,272,276,280,284,288,292,296,300,304,308/)


res@mpMinLatF     = min(tslb_xgrd@lat2d)
res@mpMaxLatF     = max(tslb_xgrd@lat2d)
res@mpMinLonF     = min(tslb_xgrd@lon2d)
res@mpMaxLonF     = max(tslb_xgrd@lon2d)
res = wrf_map_resources(f1[0],res)
res@gsnAddCyclic = False
res@tfDoNDCOverlay   = True          ; Tell NCL you are doing a native plot




plots1= new(8,graphic)
plots1(0) = gsn_csm_contour_map(wks1,daily_avg_tslb1_t(0,:,:),res)
plots1(1) = gsn_csm_contour_map(wks1,daily_avg_tslb1_t_ERA(0,:,:),res)
plots1(2) = gsn_csm_contour_map(wks1,daily_avg_tslb1_t(1,:,:),res)
plots1(3) = gsn_csm_contour_map(wks1,daily_avg_tslb1_t_ERA(1,:,:),res)
plots1(4) = gsn_csm_contour_map(wks1,daily_avg_tslb1_t(2,:,:),res)
plots1(5) = gsn_csm_contour_map(wks1,daily_avg_tslb1_t_ERA(2,:,:),res)
plots1(6) = gsn_csm_contour_map(wks1,daily_avg_tslb1_t(3,:,:),res)
plots1(7) = gsn_csm_contour_map(wks1,daily_avg_tslb1_t_ERA(3,:,:),res)

rtsP                      = True             ; modify the panel plot
rtsP@gsnMaximize          = True             ; large format
gsn_panel(wks1,plots1,((/4,2/)),rtsP)        ; draw all 'neof' as one plot


end

recite

不成熟建议：
试试看ESMF函数？

王精志

recite 发表于 2021-7-15 11:02
不成熟建议：
试试看ESMF函数？

好的，谢谢！我先看一下

王精志

recite 发表于 2021-7-15 11:02
不成熟建议：
试试看ESMF函数？

使用ESMF画出来还是这样的

王精志

recite 发表于 2021-7-15 11:02
不成熟建议：
试试看ESMF函数？

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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"

begin
path1 = "/Volumes/UNTITLED/1979/"
p1    = systemfunc("cd " + path1 +"; ls wrfout_d01_1979-0*")
;print(p1)
d = dimsizes(p1)
f1    = addfiles(path1+p1, "r")
ListSetType(f1, "cat")
printVarSummary(f1)

f = addfile ("/Volumes/UNTITLED/1979/1979.nc", "r")
f2      = addfile ("/Volumes/UNTITLED/1979/sfc.to.nc", "r")
LAT= wrf_user_getvar(f1,"XLAT",0)   ; latitude/longitude
printVarSummary(LAT)
LON = wrf_user_getvar(f1,"XLONG",0)  ; required for plotting

;lon = f->longitude({min(LON):max(LON)})
;lat = f->latitude({min(LAT):max(LAT)})

longitude = lonFlip(f->longitude)
latitude = f->latitude

lon = longitude({min(LON):max(LON)})
lat = latitude({min(LAT):max(LAT)})

;土壤温度
TSLB = wrf_user_getvar(f1,"TSLB",-1)
tslb1 = TSLB (:,0,:,:)

    Opt                   = True

;---"bilinear" is the default. "patch" and "conserve" are other options.
    Opt@InterpMethod      = "bilinear"        ;;---Change (maybe)

    Opt@WgtFileName       = "WRF_to_rect.nc"

    Opt@SrcGridLat        = LAT          ; source grid
    Opt@SrcGridLon        = LON
    Opt@SrcRegional       = True              ;;--Change (maybe)
    Opt@SrcInputFileName  = f1[0]          ; optional, but good idea

    Opt@DstGridLat        = lat           ; destination grid
    Opt@DstGridLon        = lon
    Opt@DstRegional       = True              ;;--Change (maybe)

    Opt@ForceOverwrite    = True
    Opt@PrintTimings      = True
    Opt@Debug             = True

    var_regrid = ESMF_regrid(tslb1,Opt)     ; Do the regridding

    printVarSummary(var_regrid)
   daily_avg_tslb1 = dim_avg_n_Wrap(var_regrid,(/0/))
;----------------------------------------------------------------------
; Plotting section
;
; This section creates filled contour plots of both the original
; data and the regridded data, and panels them.
;----------------------------------------------------------------------
    daily_avg_tslb1@lat2d = LAT     ; Needed for plotting. "var_regrid"
    daily_avg_tslb1@lon2d = LON     ; already has these attrs attached.

    wks = gsn_open_wks("png","/Volumes/UNTITLED/1979/spin_up_1979(TSLB)_map")

    res                       = True

    res@gsnMaximize           = True

;    res@gsnDraw               = False
;   res@gsnFrame              = False

    res@cnFillOn              = True
    res@cnLinesOn             = False
    res@cnLineLabelsOn        = False
;   res@cnFillMode            = "RasterFill"

    res@lbLabelBarOn          = False    ; Turn on later in panel

    res@mpMinLatF             = min(LAT)-1
    res@mpMaxLatF             = max(LAT)+1
    res@mpMinLonF             = min(LON)-1
    res@mpMaxLonF             = max(LON)+1

    mnmxint = nice_mnmxintvl( min(daily_avg_tslb1), max(daily_avg_tslb1), 18, False)
    res@cnLevelSelectionMode = "ManualLevels"
    res@cnMinLevelValF       = mnmxint(0)
    res@cnMaxLevelValF       = mnmxint(1)
    res@cnLevelSpacingF      = mnmxint(2)

;---Resources for plotting regridded data
    res@gsnAddCyclic  = False
    ;res@tiMainString  = "T2M regridded to a rectilinear grid (" + Opt@InterpMethod + ")"
res@tfDoNDCOverlay   = True          ; Tell NCL you are doing a native plot
    plot_regrid = gsn_csm_contour_map(wks,daily_avg_tslb1,res)


end