[ncl-talk] why the landcover index cannot display to match the colorbar?
dyjbean soybean
dyjbean at gmail.com
Wed Sep 23 20:43:30 MDT 2015
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150924/6fd61af5/attachment.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: plot_beijing_liuhuan.png
Type: image/png
Size: 100509 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150924/6fd61af5/attachment.png
More information about the ncl-talk
mailing list