meganout-月均-最大并排绘图

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"
wrffile=addfile("/data1/loucx/COPY/WRF/WRF/run/wrfout_d03_2020-08-01_00:00:00","r")
lat2d=wrffile->XLAT(0,3:176,3:170)
lon2d=wrffile->XLONG(0,3:176,3:170)
dims = dimsizes(lon2d)
geofile=addfile("/data1/loucx/COPY/CNAQ/data1/cctm/ZJ_CMAQ_appmap.nc","r")
flag=geofile->ZJ_1(0,0,:,:)
DataDir1 = "/data1/loucx/COPY/CNAQ/MEGAN/Output/"
flist1 = systemfunc("ls "+DataDir1+"megan_output_202008*_d03.ncf")
Sfile1 = addfiles(flist1, "r")
ListSetType (Sfile1, "join")
isop=Sfile1[:]->ISOPRENE
Isop=dim_avg_n(isop(:,:,0,:,:),0)
ISop=dim_avg_n(Isop(:,:,:),0)*68.12*3.6/1000  ;Dimensions and sizes:[183] x [183] t/h
do i=0,173
 do j=0,167
   if(ISop(i,j) .lt. 0) then
    ISop(i,j)=0
   end if
 end do
end do
ISop!0="row"
ISop!1="col"
ISop2=dim_avg_n(isop(:,5,0,:,:),0)*68.12*3.6/1000   ;t/h
ISop2!0="row"
ISop2!1="col"
num1=0
w=new((/29232/),float)
o=new((/29232/),float)
do z=0,173
 do f=0,167
  if(flag(z,f) .ne. 0) then
  w(num1)=ISop2(z,f) 
  o(num1)=ISop(z,f)
  num1=num1+1
  end if
  end do
end do
Min1=min(w)
Max1=max(w)
Mean1=sum(w)/num1
Min2=min(o)
Max2=max(o)
Mean2=sum(o)/num1
txid_tr = new(2,graphic)
amid_tr = new(2,graphic)
txid_tl = new(2,graphic)
amid_tl = new(2,graphic)
  txres                       = True
  txres@txPerimOn             = False ;True
  txres@txFontHeightF         = 0.028
  amres_tr                  = True
  amres_tr@amParallelPosF   =  0.48    ; This is the right edge of the plot.
  amres_tr@amOrthogonalPosF =  0.47    ; This is the top edge of the plot.
  amres_tr@amJust           = "BottomRight"
  amres_tl                  = True
  amres_tl@amParallelPosF   =  0.48    ; This is the left edge of the plot.
  amres_tl@amOrthogonalPosF = -0.48    ; This is the top edge of the plot.
  amres_tl@amJust           = "TopRight"

plot = new(2,graphic)
 wtype="png"
 wtype@wkWidth=1200
 wtype@wkHeight=1200
 wks  = gsn_open_wks (wtype,"panel2")       
 gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
 ;gsn_reverse_colormap(wks)
 res                        = True               ; plot mods desired
 res@gsnDraw                =False
 res@gsnFrame               =False
 res@gsnMaximize            = True
 res@tfDoNDCOverlay         = True ;;;;很可能是这行的原因
 res@mpDataSetName          = "Earth..4"
 res@mpDataBaseVersion = "MediumRes" 
 res@mpDataResolution = "Medium"
 res@mpProjection           = "LambertConformal"
 res@mpLambertParallel1F    = 30
 res@mpLambertParallel2F    = 60
 res@mpLambertMeridianF     = 118
 res@trGridType = "TriangularMesh"

 res@mpLimitMode            = "Corners"
 res@mpLeftCornerLatF       = lat2d(0,0)
 res@mpLeftCornerLonF       = lon2d(0,0)
 res@mpRightCornerLatF      = lat2d(dims(0)-1,dims(1)-1)
 res@mpRightCornerLonF      = lon2d(dims(0)-1,dims(1)-1)
 res@mpGeophysicalLineThicknessF = 0.6 ;设置用于绘制地球物理边界轮廓的线的厚度
 res@pmTickMarkDisplayMode = "Always"
 res@tmXTOn = False ;不要绘制顶部轴刻度线
 res@tmYROn = False ;不要绘制右边轴刻度线
 res@tmXBLabelFontHeightF = 0.02
 res@cnFillOn = True
 res@cnLinesOn = False
 res@cnLevelSelectionMode = "ExplicitLevels"
 res@cnLineLabelsOn         = False
 res@lbLabelBarOn = True  ;开启图形各自的色标
 res@lbTopMarginF = 0.15
 res@pmLabelBarHeightF = 0.09
 res@pmLabelBarOrthogonalPosF = 0.02  ;上下移动
 res@lbLabelFontHeightF  = 0.0225
 res@lbTitleFontHeightF = 0.025
 res@lbTitleOffsetF = 0.12
 res@lbBoxLinesOn = True
 res@tiMainString           = "Average emission rate in August"
 res@tiMainFontHeightF      = 0.035              ; smaller title
 res@tiMainOffsetYF         = -0.005
 res@gsnLeftString              = "t/h"             
 res@gsnLeftStringFontHeightF=0.035
 res@gsnRightString             = ""
;  res@mpFillOn             = True
;  res@mpOceanFillColor     = "White"
;  res@mpLandFillColor      = "transparent"
;  res@mpFillDrawOrder      = "PostDraw"
 plres = True
 plres@gsnMaximize       = True
 plres@amJust   = "TopLeft"
 plres@gsnFrame         = False
 plres@gsnPanelDebug = False
 plres@gsnPanelBoxes = False
 plres@gsnPanelFigureStringsPerimOn = False
 plres@gsnPanelFigureStringsFontHeightF = 0.013
 plres@gsnPanelFigureStringsBackgroundFillColor = "transparent"


 res1                        = True 
 res1@cnLevels=fspan(0,0.043,44)
 res1@gsnSpreadColors=True
 res1@gsnSpreadColorStart=0
 res1@lbLabelStride = 10
 res1@gsnCenterString  ="a"

 res2                        = True
 res2@cnLevels=fspan(0,0.15,16)
 res2@gsnSpreadColors=True
 res2@gsnSpreadColorStart=0
 res2@lbLabelStride = 2
 res2@gsnCenterString  ="b"


 shp1="/data1/loucx/Build_WRF/cnmap/浙江/Zhejiang.shp"
 shp2="/data1/loucx/Build_WRF/cnmap/浙江/CHN_adm1.shp"

 dumstr1 = unique_string("marker1")
 plot(0) = gsn_csm_contour_map(wks,ISop,res)
  ;tl_label = "a " ;;;;;;;
  ;txid_tl(0) = gsn_create_text(wks, tl_label, txres)
  ;amid_tl(0) = gsn_add_annotation(plot(0), txid_tl(0), amres_tl)
  tr_label = "max=" + sprintf("%5.2f",Max2) + " min=" + sprintf("%5.2f",Min2)+ " mean=" + sprintf("%5.2f",Mean2)
 txid_tr(0) = gsn_create_text(wks,tr_label, txres)
 amid_tr(0) = gsn_add_annotation(plot(0), txid_tr(0), amres_tr)
 lnres1        = True     
 lnres1@gsLineColor      ="Black"
 lnres1@gsLineThicknessF = 1.2
 shp_plot0 = gsn_add_shapefile_polylines(wks,plot(0),shp1,lnres1)
 lnres2        = True
 lnres2@gsLineColor      ="Black"
 lnres2@gsLineThicknessF = 2         ; 2x thickness  ---2 to 3.5
 shp_plot1 = gsn_add_shapefile_polylines(wks,plot(0),shp2,lnres2)
;-------------------------------
res@tiMainString           = "13: 00 emission rate"

 dumstr2 = unique_string("marker1")
 plot(1) = gsn_csm_contour_map(wks,ISop2,res)
  ;tl_label = "b " ;;;;;;;
  ;txid_tl(1) = gsn_create_text(wks, tl_label, txres)
  ;amid_tl(1) = gsn_add_annotation(plot(1), txid_tl(1), amres_tl)
 tr_label = "max=" + sprintf("%5.2f",Max1) + " min=" + sprintf("%5.2f",Min1)+ " mean=" + sprintf("%5.2f",Mean1)
 txid_tr(1) = gsn_create_text(wks,tr_label, txres)
 amid_tr(1) = gsn_add_annotation(plot(1), txid_tr(1), amres_tr)

 lnres3        = True
 lnres3@gsLineColor      ="Black"
 lnres3@gsLineThicknessF = 1.2
 shp_plot2 = gsn_add_shapefile_polylines(wks,plot(1),shp1,lnres3)
 lnres4        = True
 lnres4@gsLineColor      ="Black"
 lnres4@gsLineThicknessF = 2         ; 2x thickness  ---2 to 3.5
 shp_plot3 = gsn_add_shapefile_polylines(wks,plot(1),shp2,lnres4)
;;;去白边
 boxlat=(/lat2d(0,165),lat2d(0,167),lat2d(173,167),lat2d(173,165),lat2d(0,165)/)
 boxlon=(/lon2d(0,165),lon2d(0,167),lon2d(173,167),lon2d(173,165),lon2d(0,165)/)
 gonres=True
 gonres@gsFillColor="White"
 dum1=gsn_add_polygon(wks,plot(0),boxlon,boxlat,gonres)
 dum2=gsn_add_polygon(wks,plot(1),boxlon,boxlat,gonres);;;添加矩形去掉边 v

 pres                        = True
 pres@gsnMaximize            = True
 gsn_panel(wks,plot,(/1,2/),pres)