[ncl-talk] why the landcover index cannot display to match the colorbar?
Mary Haley
haley at ucar.edu
Thu Sep 24 09:50:25 MDT 2015
First: you want to make sure gsnSpreadColors is *not* set at all. When you
use cnFillPalette, you shouldn't need the gsnSpreadXXXX resources.You've
already made sure you have the exact number of colors that you need for
your contour levels, so you're good there.
The issue is that NCL always contours the data based on values that fall
between the given levels, and not exactly at the given levels. So, for
example, if you have data values that fall exactly at values 1, 2, 3, 4, 5,
then to create contours at these values, you need to use values that fall
between them:
res at cnLevels = (/1.5, 2.5, 3.5, 4.5 /)
With four contour levels, you will get five colors in your labelbar:
- the first color in your labelbar represents all values <1.5 (hence
covering the value 1)
- the second color will represent all values >= 1.5 and < 2.5 (hence
covering the value 2)
. . .
- the last (fifth) color will represent all values > 4.5
If you look at the original vegland_1.ncl script, it is setting the levels
to go from 2 to 18 in steps of 1, but the actual label strings go from 1 to
18 in steps of one. This works because then the first color then
represents all values < 2 (thus catching the values equal to 1), and so on.
I just like to use values like 1.5, 2.5, etc, to make it more clear, but
the end result is the same.
I think the rest of your resources for setting up centered-labels and
setting the label strings themselves looks good.
One final comment: gsn_panel only rearranges the existing plots and adds a
labelbar, if requested. It does *not* change the existing plots. Thus, the
"cnLevel" resource will not be recognized by gsn_panel.
Also, to make sure you get every box labeled, make sure you set:
pres at lbLabelStride = 1
pres at lbLabelAutoStride = False
--Mary
On Wed, Sep 23, 2015 at 8:43 PM, dyjbean soybean <dyjbean at gmail.com> wrote:
> hi,
> i have made a landcover map, which includes indices from 1 to 19 and
> 31-32 , all is 22 categories.
>
> but after i ploted them, only the 1 to 19 and 33 category appearing in
> figure, but the 31-32 these two category cannot be shown in figure.
>
> those below is my NCL script:
>
>
> ++++++++++++++++++++++++++++++++++
>
> 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="globecover_hubei_liuhuan_usgs_ronghe_snow-ice1.txt"
> fconname="globecover_huabei_liuhuan_usgs_new_allocate1.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 at _FillValue)
> delete(ndat at _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 at units="degrees_north"
> olat at long_name="latitude"
>
> olon!0="lon"
> olon at units="degrees_east"
> olon at long_name="longitude"
>
> odat!0="lat"
> odat!1="lon"
> odat&lat=olat
> odat&lon=olon
> odat at units=""
> odat at long_name="land cover confused with globaland30 and globercover"
> odat at _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 at units="degrees_north"
> nlat at long_name="latitude"
>
> nlon!0="lon"
> nlon at units="degrees_east"
> nlon at long_name="longitude"
>
> ndat!0="lat"
> ndat!1="lon"
> ndat&lat=nlat
> ndat&lon=nlon
> ndat at units=""
> ndat at long_name="land cover new categories after confused"
> ndat at _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 at return_mask = True
> ; opto at debug = True
> opto at minlon = minlon ; Makes the shapefile masking
> opto at maxlon = maxlon ; go faster.
> opto at minlat = minlat
> opto at 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 at return_mask = True
> ; optn at debug = True
> optn at minlon = minlon ; Makes the shapefile masking
> optn at maxlon = maxlon ; go faster.
> optn at minlat = minlat
> optn at 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 at _FillValue)
> copy_VarMeta(odat,omask)
>
> ;---Create masked data array
> nmask = where(ndat_mask.eq.1,ndat,ndat at _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 at gsnDraw = False
> res at gsnFrame = False
> ;res at gsnMaximize = True ; ps, pdf
> res at gsnAddCyclic = False
>
> ;res at mpMinLatF = min(odat&lat)
> ;res at mpMaxLatF = max(odat&lat)
> ;res at mpMinLonF = min(odat&lon)
> ;res at mpMaxLonF = max(odat&lon)
>
> res at cnFillOn = True ; color Fill
> res at cnFillMode = "RasterFill"
> res at cnLinesOn = False ; Turn off contour lines
> res at cnLineLabelsOn = False
> res at cnLevelSelectionMode = "ExplicitLevels" ; set explict contour levels
> res at cnLevels = tobyte(ispan(1,ninfo-1,1)+1)
> res at cnFillPalette = colorscheme ; distinct colors for
> categories
> res at gsnSpreadColors = False ; use each color sequentially
> res at cnInfoLabelOn = False
>
>
> res at mpFillOn = False
>
>
> res at gsnAddCyclic = False ; regional data
> ;res at trGridType = "TriangularMesh"
> res at tmXBLabelFontHeightF = 0.01 ; Make lon & lat text smaller
> res at tmYLLabelFontHeightF = res at tmXBLabelFontHeightF
>
> res at gsnLeftString = ""
> res at gsnRightString = ""
>
> res1 = True
> res1 = res
> res1 at mpMinLatF = minlat
> res1 at mpMaxLatF = maxlat
> res1 at mpMinLonF = minlon
> res1 at mpMaxLonF = maxlon
>
> res1 at tmXBMode = "Explicit"
> res1 at tmXBValues =
> (/116.0,116.1,116.2,116.3,116.4,116.5,116.6,116.7,116.8/)
> res1 at 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 at tmYLMode = "Explicit"
> res1 at tmYLValues = (/39.6,39.7,39.8,39.9,40.0,40.1,40.2/)
> res1 at 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 at tiMainString = "confused LC"
> res1 at lbLabelBarOn = False
> plot(0) = gsn_csm_contour_map_ce(wks, omask, res1) ; create plot
>
> res2 = True
> res2=res
> res2 at mpMinLatF = minlat
> res2 at mpMaxLatF = maxlat
> res2 at mpMinLonF = minlon
> res2 at mpMaxLonF = maxlon
> res2 at lbLabelBarOn = False
>
> res2 at tmXBMode = "Explicit"
> res2 at tmXBValues =
> (/116.0,116.1,116.2,116.3,116.4,116.5,116.6,116.7,116.8/)
> res2 at 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 at tmYLMode = "Explicit"
> res2 at tmYLValues = (/39.6,39.7,39.8,39.9,40.0,40.1,40.2/)
> res2 at 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 at tiMainString = "reclassified LC"
> plot(1) = gsn_csm_contour_map_ce(wks,nmask,res2)
>
> pres = True
> pres at gsnFrame = False
> pres at gsnPanelLabelBar = True
> pres at lbLabelFontHeightF = 0.01
> pres at gsnPanelBottom = 0.05
> pres at lbLabelPosition = "Center" ; label position
> pres at lbLabelAlignment = "BoxCenters" ; label orientation
> pres at cnLevels = (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,\
> 31,32,33/)
> pres at lbLabelStrings = (/"1","2","3","4","5","6","7","8","9","10",\
> "11","12","13","14","15","16","17","18","19",\
> "31","32","33"/)
> ;pres at lbLabelStrings = ispan(0,ninfo,1)
> ;pres at lbLabelStride = 1
> pres at lbLabelAutoStride = True
>
> gsn_panel(wks,plot,(/1,2/),pres)
>
> rtxt = True
> rtxt at txJust = "TopLeft"
> rtxt at 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 at txFontHeightF
> end if
> end do
> xx = xx + 0.20
> end do
>
> ;draw(plot)
> frame(wks)
>
> end
>
>
>
> ++++++++++++++++++++++++++++++++++
>
> i don't know if the nested pic can show, so i attach it
>
> i change the urban (denote with red color in left panel) into impervious
> category (31,32,33，show in right panel), but eventually only 33 appeared,
> 31 and 32 disappeared.
> the 31-32 exists in data.
>
> can the ncl script match the landcover index automatically?
>
>
> would you please give me some advices for this problem which i have been
> confused some days？
>
> thanks
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150924/27ca6f5d/attachment.html
More information about the ncl-talk
mailing list