在画landcover时遇到问题，colorbar指定的颜色在图中没有显示

novelbean

hi，
  你好！

  在画经过重新分类landcover数据时，数据中有31-33新的分类，在colorbar中也有指定颜色，但是在图中没有全部显示，只显示了33这种分类的颜色，请问是什么原因？
  以下是脚本和图
  ++++++++++++++++++++++++++++++++++++++++
  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/csm/shea_util.ncl"
load "./shapefile_mask_data.ncl"
begin
forgname="old.txt"
fconname="new.txt"
onlat=5403
omlon=5403
nnlat=1801
nmlon=1801
odat=new((/onlat,omlon/),"integer","No_FillValue")
ndat=new((/nnlat,nmlon/),"integer","No_FillValue")
odat=readAsciiTable(forgname,omlon,"integer",5)
ndat=asciiread(fconname,(/nnlat,nmlon/),"integer")
delete(odat@_FillValue)
delete(ndat@_FillValue)
;; confused lc
olat=new((/onlat/),"float","No_FillValue")
olon=new((/omlon/),"float","No_FillValue")
do i=0,onlat-1
olat(i)=45.00413067-i*0.0027778
end do

do j=0,omlon-1
olon(j)=109.99606533+j*0.0027778
end do
olat!0="lat"
olat@units="degrees_north"
olat@long_name="latitude"
olon!0="lon"
olon@units="degrees_east"
olon@long_name="longitude"
odat!0="lat"
odat!1="lon"
odat&lat=olat
odat&lon=olon
odat@units=""
odat@long_name="land cover confused with globaland30 and globercover"
odat@_FillValue = -999
;; new classified lc
nlat=new((/nnlat/),"float","No_FillValue")
nlon=new((/nmlon/),"float","No_FillValue")
do i=0,nnlat-1
nlat(i)=45.00413067-i*0.0083334
end do

do j=0,nmlon-1
nlon(j)=109.99606533+j*0.0083334
end do
nlat!0="lat"
nlat@units="degrees_north"
nlat@long_name="latitude"
nlon!0="lon"
nlon@units="degrees_east"
nlon@long_name="longitude"
ndat!0="lat"
ndat!1="lon"
ndat&lat=nlat
ndat&lon=nlon
ndat@units=""
ndat@long_name="land cover new categories after confused"
ndat@_FillValue = -999
print(odat(onlat-1,omlon-1))
;
; If you already have the mask NetCDF file, set this to False.
; Creating the mask can be slow!
;
   CREATE_MASK = True
;---Whether to draw lat/lon points and shapefile outlines
   ADD_LATLON_POINTS      = False
   ADD_SHAPEFILE_OUTLINES = True
;---Name of shapefile containing USA outlines
   shp_fname = "huan_6.shp"
;---Name of file to write mask to or to read mask from.
   mask_oname = "liuhuan_mask_300m.nc"
   mask_nname = "liuhuan_mask_1km.nc"
;---Rough area we are interested in. Everything outside this will be masked.
   minlon = 116.0
   maxlon = 116.8
   minlat =   39.6
   maxlat =   40.2
;; odat
if(CREATE_MASK) then
     print("Creating the mask file...")
;---Create a new mask using a shapefile of USA
     odims = dimsizes(odat)
     opto             = True
     opto@return_mask = True
;    opto@debug       = True
     opto@minlon      = minlon     ; Makes the shapefile masking
     opto@maxlon      = maxlon     ; go faster.
     opto@minlat      = minlat
     opto@maxlat      = maxlat
     odat_mask        = shapefile_mask_data(odat,shp_fname,opto)
;---Write new mask to file
     system("rm -f " + mask_oname)
     oout           = addfile(mask_oname,"c")
     oout->odat_mask = odat_mask
   else
     print("Reading mask off file.")
;---Read the new mask from the NetCDF file
     omask    = addfile(mask_oname,"r")
     odat_mask = omask->odat_mask
   end if
;; ndat
if(CREATE_MASK) then
     print("Creating the mask file...")
;---Create a new mask using a shapefile of USA
     ndims = dimsizes(ndat)
     optn             = True
     optn@return_mask = True
;    optn@debug       = True
     optn@minlon      = minlon     ; Makes the shapefile masking
     optn@maxlon      = maxlon     ; go faster.
     optn@minlat      = minlat
     optn@maxlat      = maxlat
     ndat_mask        = shapefile_mask_data(ndat,shp_fname,optn)
;---Write new mask to file
     system("rm -f " + mask_nname)
     nout           = addfile(mask_nname,"c")
     nout->ndat_mask = ndat_mask
   else
     print("Reading mask off file.")
;---Read the new mask from the NetCDF file
     nmask    = addfile(mask_nname,"r")
     ndat_mask = nmask->ndat_mask
   end if
;---Create masked data array
   omask = where(odat_mask.eq.1,odat,odat@_FillValue)
   copy_VarMeta(odat,omask)
;---Create masked data array
   nmask = where(ndat_mask.eq.1,ndat,ndat@_FillValue)
   copy_VarMeta(ndat,nmask)

  ;ind32=ind(ndtooned(nmask).eq.32)
  ;print(ind32)
  dimsnmk=dimsizes(nmask)
  do i=0,dimsnmk(0)-1
    do j=0,dimsnmk(1)-1
       if(any(ndat(i,j).eq.32)) then
          print(i+" "+j+" "+" "+ndat(i,j))
       end if
    end do
  end do
;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; create plot
;;;;;;;;;;;;;;;;;;;;;;;;;;;;
plot=new(2,"graphic")
info=(/                        \
      "1 urban_and_built-in"   ,\
      "2 dryland_cropland"     ,\
      "3 irrigated_cropland"   ,\
      "4 mixed_dry_cropland"   ,\
      "5 cropland_grassland"   ,\
      "6 cropland_woodland"    ,\
      "7 grassland"            ,\
      "8 shrubland"            ,\
      "9 mixed_shrub_grassland",\
      "10 savanna"              ,\
      "11 deciduous_broadleaf"  ,\
      "12 deciduous_needleleaf" ,\
      "13 evergreen_broadleaf"  ,\
      "14 evergreen_needleleaf" ,\
      "15 mixed_forest"         ,\
      "16 water"                ,\
      "17 herbaceous_wetland"   ,\
      "18 wooded_wetland"       ,\
      "19 barren_or_sparse"     ,\
      "31 low_intensity_res"    ,\
      "32 mid_intensity_res"    ,\
      "33 high_intensity_res"/)
ninfo=dimsizes(info)
  ninfo = dimsizes(info)        ; # of classifications

;************************************************
; create plot
;************************************************
  pltType = "png"               ; ps, pdf, png, x11, eps

  colorscheme =(/"red","lawngreen","mediumseagreen","springgreen",\
       "limegreen","darkolivegreen3","green","lightsalmon","green3",\
        "gold","palegreen","olivedrab3","chartreuse3","chartreuse",       \
        "darkgreen","blue","hotpink2","hotpink","light cyan",           \
        "khaki1","violetred1","purple3"/)
  ncolors = dimsizes(colorscheme)
  if (ninfo.ne.ncolors) then             ; make sure # of colors match categories (classes)
      print("size mismatch: ninfo="+ninfo+"   ncolors="+ncolors)
      exit
  end if
  wks = gsn_open_wks(pltType,"plot_origin_new_lc_test")
  plot=new(2,graphic)
  res                  = True                ; plot mods desired
  res@gsnDraw          = False
  res@gsnFrame         = False
  ;res@gsnMaximize      = True               ; ps, pdf
  res@gsnAddCyclic     = False
  ;res@mpMinLatF   = min(odat&lat)
  ;res@mpMaxLatF   = max(odat&lat)
  ;res@mpMinLonF   = min(odat&lon)
  ;res@mpMaxLonF   = max(odat&lon)
  res@cnFillOn         = True               ; color Fill
  res@cnFillMode       = "RasterFill"
  res@cnLinesOn        =  False             ; Turn off contour lines
  res@cnLineLabelsOn   =  False
  res@cnLevelSelectionMode = "ExplicitLevels" ; set explict contour levels
  res@cnLevels         = tobyte(ispan(1,ninfo-1,1)+1)
  res@cnFillPalette    = colorscheme        ; distinct colors for categories
  res@gsnSpreadColors  = False              ; use each color sequentially
  res@cnInfoLabelOn   = False

  res@mpFillOn         = False

  res@gsnAddCyclic   = False             ; regional data
  ;res@trGridType     = "TriangularMesh"
  res@tmXBLabelFontHeightF = 0.01        ; Make lon & lat text smaller
  res@tmYLLabelFontHeightF = res@tmXBLabelFontHeightF
  res@gsnLeftString  = ""
  res@gsnRightString = ""
  res1  =  True
  res1  =  res
  res1@mpMinLatF   = minlat
  res1@mpMaxLatF   = maxlat
  res1@mpMinLonF   = minlon
  res1@mpMaxLonF   = maxlon
  res1@tmXBMode    = "Explicit"
  res1@tmXBValues   = (/116.0,116.1,116.2,116.3,116.4,116.5,116.6,116.7,116.8/)
  res1@tmXBLabels  =  (/"116.0~S~o~N~E","116.1~S~o~N~E","116.2~S~o~N~E","116.3~S~o~N~E","116.4~S~o~N~E","116.5.0~S~o~N~E",\
                        "116.6~S~o~N~E","116.7~S~o~N~E","116.8~S~o~N~E"/)
  res1@tmYLMode    = "Explicit"
  res1@tmYLValues   = (/39.6,39.7,39.8,39.9,40.0,40.1,40.2/)
  res1@tmYLLabels  = (/"39.6~S~o~N~N","39.7~S~o~N~N","39.8~S~o~N~N","39.9~S~o~N~N","40.0~S~o~N~N","40.1~S~o~N~N","41.2~S~o~N~N"/)
  res1@tiMainString   = "confused LC"
  res1@lbLabelBarOn   =  False
  plot(0) = gsn_csm_contour_map_ce(wks, omask, res1) ; create plot
  res2 = True
  res2=res
  res2@mpMinLatF   = minlat
  res2@mpMaxLatF   = maxlat
  res2@mpMinLonF   = minlon
  res2@mpMaxLonF   = maxlon
  res2@lbLabelBarOn   =  False
  res2@tmXBMode    = "Explicit"
  res2@tmXBValues  = (/116.0,116.1,116.2,116.3,116.4,116.5,116.6,116.7,116.8/)
  res2@tmXBLabels  = (/"116.0~S~o~N~E","116.1~S~o~N~E","116.2~S~o~N~E","116.3~S~o~N~E","116.4~S~o~N~E","116.5.0~S~o~N~E",\
                       "116.6~S~o~N~E","116.7~S~o~N~E","116.8~S~o~N~E"/)
  res2@tmYLMode    = "Explicit"
  res2@tmYLValues   = (/39.6,39.7,39.8,39.9,40.0,40.1,40.2/)
  res2@tmYLLabels  = (/"39.6~S~o~N~N","39.7~S~o~N~N","39.8~S~o~N~N","39.9~S~o~N~N","40.0~S~o~N~N","40.1~S~o~N~N","41.2~S~o~N~N"/)
  res2@tiMainString =  "reclassified LC"
  plot(1) =  gsn_csm_contour_map_ce(wks,nmask,res2)
  pres = True
  pres@gsnFrame =  False
  pres@gsnPanelLabelBar  = True
  pres@lbLabelFontHeightF = 0.01
  pres@gsnPanelBottom = 0.05
  pres@lbLabelPosition  = "Center"           ; label position
  pres@lbLabelAlignment = "BoxCenters"       ; label orientation
  pres@cnLevels   =  (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,\
                       31,32,33/)
  pres@lbLabelStrings    = (/"1","2","3","4","5","6","7","8","9","10",\
                           "11","12","13","14","15","16","17","18","19",\
                           "31","32","33"/)
;pres@lbLabelStrings   = ispan(0,ninfo,1)
  ;pres@lbLabelStride    = 1
  pres@lbLabelAutoStride = True
  gsn_panel(wks,plot,(/1,2/),pres)
  rtxt = True
  rtxt@txJust        = "TopLeft"
  rtxt@txFontHeightF = 0.0125
  ; Add text: rows x columns of text (arbitrary)
; Usually must play with xx, yy and txFontHeightF
  nrow = 5       ; # rows
  ncol = 6       ; # columns
  n  = -1
  xx = 0.015                ; arbitrary
  do nc=0,ncol-1
     yy = 0.20
    do nr=0,nrow-1
       n = n+1
       if (n.le.(ninfo-1)) then
           gsn_text_ndc (wks,info(n),xx,yy,rtxt)
           yy = yy - 2*rtxt@txFontHeightF
       end if
    end do
    xx = xx + 0.20
  end do
;draw(plot)
frame(wks)
end

  ++++++++++++++++++++++++++++++++++++++++

问题是为什么数据中的31-33这三种分类，在图中体现不出来，只显示了一种colorbar中给定的33这种分类

请指教，谢谢

novelbean

知道原因了，因为分类不是连续的，有离散的分类
所以
res@cnLevels         = tobyte((/2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,31,32,33/)) ;tobyte(ispan(1,ninfo-1,1)+1)

要改成这样子，才能正确显示

novelbean

ncl中cnLevels属性应该根据赋的值来与数据中像匹配，然后再做处理

Genre

{:5_213:}{:5_213:}{:5_213:}

lgz

{:5_213:}