如何用NCL提出NC文件某个经纬度的温度或者NO浓度信息与观测点的作比较

luo981858

这是一个画出一个CMAQ排除的NC文件，模拟的一个区域每个网格NO的浓度分布。
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 "/data4/huihong/map_info/cnmap/cnmap.ncl"

begin

  ;grid file
  grid_fn="/data4/huihong/map_info/cmaq_domain/GRIDCRO2D_3km"

  ;Fig title
  fig_title="NO Concentration (ppmv)"

  g_font=0.02

  ;File IO
  pdata_fn     = "/data4/huihong/data/cctm/multi-nodes/CCTM_e1a_Linux2_x86_64intel.ACONC.CN3GD_152X110.multi-nodes.2017232"
  ctrl_in = addfile(pdata_fn ,"r")  ; open output netCDF file
  var1    = ctrl_in->NO(:,0,:,:)
; var1    = var1+ctrl_in->NO2(:,0,:,:)
  printVarSummary(var1)
  print(max(var1))
  print(min(var1))
  latlon_in   =  addfile(grid_fn,"r")
  lat = latlon_in->LAT(0,0,:,:)
  lon = latlon_in->LON(0,0,:,:)

  var1@lat2d=lat
  var1@lon2d=lon

  res                  = True       ; use plot options

  res@cnFillOn             = True               ; turn on color for contours
  res@cnLinesOn            = False              ; turn off contour lines
;  res@cnFillOn             = True               ; turn on color for contours
  res@cnLinesOn            = False              ; turn off contour lines
  res@cnLineLabelsOn       = False              ; turn off contour line labels

;  res@cnLevelSelectionMode = "ExplicitLevels"   ; set manual contour levels
;  res@cnLevels = ispan(0,60,10)
;  res@cnFillColors =(/-1,20,47,57,94,127,152/)

res@cnLevelSelectionMode  = "ManualLevels" ; set manual contour levels
res@cnMinLevelValF        = 0          ; set min contour level
res@cnMaxLevelValF        = 0.06          ; set max contour level
res@cnLevelSpacingF       = 0.002          ; set contour interval

  res@gsnFrame         = False
  res@gsnDraw  = False

  res@gsnSpreadColors      = True               ; use full color map
;  res@gsnSpreadColorStart  = 2               ; start at color 17
;  res@gsnSpreadColorEnd    = 14                ; end at color 200

  res@gsnRightString = ""
  res@gsnMaximize      = True       ; fill up the page
  ;res@gsnAddCyclic   = True;False
  res@gsnPaperOrientation = "portrait"
  res@gsnContourZeroLineThicknessF = 2.  ;set thickness of zero

  res@cnFillMode           = "CellFill" ; Raster Mode

  res@lbLabelBarOn = True   ; turn off the label bar
  res@lbOrientation          = "vertical"
  res@lbLabelFontHeightF  = 0.02              ; make labels smaller
  ;res@lbLabelStride = 1

  res@mpMinLatF            = min(var1@lat2d)        ; zoom in on map

  res@mpMinLatF            = min(var1@lat2d)        ; zoom in on map
  res@mpMaxLatF            = max(var1@lat2d)
  res@mpMinLonF            = min(var1@lon2d)
  res@mpMaxLonF            = max(var1@lon2d)

  res@mpGeophysicalLineThicknessF = 2.0 ;costal line thick
  res@tmXBTickSpacingF = 2
  res@tmYLTickSpacingF = 2

  res@tmXBLabelFontHeightF =g_font
  res@tmYLLabelFontHeightF = g_font
  res@gsnStringFontHeightF = g_font
  res@tiMainFontHeightF= g_font
  ;res@lbLabelFontHeightF = 0.02
  ;res@pmLabelBarOrthogonalPosF = .12           ; move label bar down
  res@tmXBMajorThicknessF = 2.0
  res@tmYLMajorThicknessF = 2.0
  res@tmXBMinorThicknessF = 2.0
  res@tmYLMinorThicknessF = 2.0
  res@tmBorderThicknessF = 2.0
  res@tmYLMajorLengthF = 0.002


  res@mpFillOn                = True
  res@mpOutlineOn             = False  ; Use outlines from shapefile
  res@cnFillDrawOrder         = "PreDraw"
  res@mpDataBaseVersion       = "MediumRes"
  res@mpDataSetName           = "Earth..4"
  res@mpAreaMaskingOn         = True
  res@mpMaskAreaSpecifiers    = (/"China","Taiwan","Disputed area between India and China","India:Arunachal Pradesh"/)
  res@mpLandFillColor         = "white"
  res@mpInlandWaterFillColor  = "white"
  ;res@mpOceanFillColor        = "white"
  res@mpInlandWaterFillColor  = "white"
  ;res@mpOceanFillColor        = "white"
  ;res@mpOutlineBoundarySets   = "NoBoundaries"

;>============================================================<
;                      add China map
;>------------------------------------------------------------<
  cnres           = True
  cnres@china     = False       ;draw china map or not
  cnres@river     = False       ;draw changjiang&huanghe or not
  cnres@province  = True       ;draw province boundary or notcnres@nanhai    = False       ;draw nanhai or not
  cnres@nanhai    = False       ;draw nanhai or not
  cnres@diqu      = True       ; draw diqujie or not

  do i_hour=0,23
        ;hours=(/0920,0921,0922,0923/)
        g_fig_name="/data4/huihong/fig/NO"+i_hour
        print(i_hour)

        res@gsnCenterString="NO_hourly"
        res@gsnLeftString="20150820"+i_hour
        res@gsnRightString="ppmv"
        wks = gsn_open_wks("png",g_fig_name)       ; open file to plot
        ;gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
        gsn_define_colormap(wks,"wh-bl-gr-ye-re")
    ;    gsn_reverse_colormap(wks)
        plot = gsn_csm_contour_map(wks,var1(i_hour,:,:),res) ; dim1 = 0 for lvl = 0
        chinamap = add_china_map(wks,plot,cnres)
        draw(plot)
        frame(wks)

        delete([/wks,plot,chinamap/])
  end do

彭雪娇

很棒的脚本！！楼主厉害！！学到了，谢谢！！

2160409394

厉害厉害