[ncl-talk] why the landcover index cannot display to match the colorbar?
dyjbean soybean
dyjbean at gmail.com
Thu Sep 24 19:57:24 MDT 2015
thanks for your advices,following your suggestion, i changed the
res at cnLevels=res at cnLevels =
tobyte((/2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,31,32,33/))
i think the index in this attribute is corresponding to the number in data.
2015-09-24 23:50 GMT+08:00 Mary Haley <haley at ucar.edu>:
>
> 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
>>
>>
>
--
通讯地址：北京市朝阳区大屯路9718信箱
中国科学院遥感应用研究所
邮编：100101
电话：010-64860422
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150925/55876d73/attachment.html
More information about the ncl-talk
mailing list