<div dir="ltr">hi,<div> i have made a landcover map, which includes indices from 1 to 19 and 31-32 , all is 22 categories.<div><br></div><div>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.</div><div><br></div>those below is my NCL script:
</div><div><br></div><div><br></div><div>++++++++++++++++++++++++++++++++++</div><div><div><br></div><div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"</div><div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"</div><div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"</div><div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"</div><div>load "./shapefile_mask_data.ncl"</div><div><br></div><div>begin</div><div>forgname="globecover_hubei_liuhuan_usgs_ronghe_snow-ice1.txt"</div><div>fconname="globecover_huabei_liuhuan_usgs_new_allocate1.txt"</div><div><br></div><div>onlat=5403</div><div>omlon=5403</div><div><br></div><div>nnlat=1801</div><div>nmlon=1801</div><div><br></div><div>odat=new((/onlat,omlon/),"integer","No_FillValue")</div><div>ndat=new((/nnlat,nmlon/),"integer","No_FillValue")</div><div><br></div><div>odat=readAsciiTable(forgname,omlon,"integer",5)</div><div>ndat=asciiread(fconname,(/nnlat,nmlon/),"integer")</div><div>delete(odat@_FillValue)</div><div>delete(ndat@_FillValue)</div><div><br></div><div>;; confused lc</div><div>olat=new((/onlat/),"float","No_FillValue")</div><div>olon=new((/omlon/),"float","No_FillValue")</div><div><br></div><div>do i=0,onlat-1</div><div> olat(i)=45.00413067-i*0.0027778</div><div>end do</div><div><br></div><div><br></div><div>do j=0,omlon-1</div><div> olon(j)=109.99606533+j*0.0027778</div><div>end do</div><div><br></div><div>olat!0="lat"</div><div>olat@units="degrees_north"</div><div>olat@long_name="latitude"</div><div><br></div><div>olon!0="lon"</div><div>olon@units="degrees_east"</div><div>olon@long_name="longitude"</div><div><br></div><div>odat!0="lat"</div><div>odat!1="lon"</div><div>odat&lat=olat</div><div>odat&lon=olon</div><div>odat@units=""</div><div>odat@long_name="land cover confused with globaland30 and globercover"</div><div>odat@_FillValue = -999</div><div><br></div><div>;; new classified lc</div><div>nlat=new((/nnlat/),"float","No_FillValue")</div><div>nlon=new((/nmlon/),"float","No_FillValue")</div><div><br></div><div>do i=0,nnlat-1</div><div> nlat(i)=45.00413067-i*0.0083334</div><div>end do</div><div><br></div><div><br></div><div>do j=0,nmlon-1</div><div> nlon(j)=109.99606533+j*0.0083334</div><div>end do</div><div><br></div><div>nlat!0="lat"</div><div>nlat@units="degrees_north"</div><div>nlat@long_name="latitude"</div><div><br></div><div>nlon!0="lon"</div><div>nlon@units="degrees_east"</div><div>nlon@long_name="longitude"</div><div><br></div><div>ndat!0="lat"</div><div>ndat!1="lon"</div><div>ndat&lat=nlat</div><div>ndat&lon=nlon</div><div>ndat@units=""</div><div>ndat@long_name="land cover new categories after confused"</div><div>ndat@_FillValue = -999</div><div><br></div><div>print(odat(onlat-1,omlon-1))</div><div><br></div><div>;</div><div>; If you already have the mask NetCDF file, set this to False.</div><div>; Creating the mask can be slow!</div><div>;</div><div> CREATE_MASK = True</div><div><br></div><div>;---Whether to draw lat/lon points and shapefile outlines</div><div> ADD_LATLON_POINTS = False</div><div> ADD_SHAPEFILE_OUTLINES = True</div><div><br></div><div>;---Name of shapefile containing USA outlines</div><div> shp_fname = "huan_6.shp"</div><div><br></div><div>;---Name of file to write mask to or to read mask from.</div><div> mask_oname = "<a href="http://liuhuan_mask_300m.nc">liuhuan_mask_300m.nc</a>"</div><div> mask_nname = "<a href="http://liuhuan_mask_1km.nc">liuhuan_mask_1km.nc</a>"</div><div><br></div><div>;---Rough area we are interested in. Everything outside this will be masked.</div><div> minlon = 116.0</div><div> maxlon = 116.8</div><div> minlat = 39.6</div><div> maxlat = 40.2</div><div><br></div><div>;; odat</div><div>if(CREATE_MASK) then</div><div> print("Creating the mask file...")</div><div><br></div><div>;---Create a new mask using a shapefile of USA</div><div> odims = dimsizes(odat)</div><div> opto = True</div><div> opto@return_mask = True</div><div>; opto@debug = True</div><div> opto@minlon = minlon ; Makes the shapefile masking</div><div> opto@maxlon = maxlon ; go faster.</div><div> opto@minlat = minlat</div><div> opto@maxlat = maxlat</div><div> odat_mask = shapefile_mask_data(odat,shp_fname,opto)</div><div><br></div><div>;---Write new mask to file</div><div> system("rm -f " + mask_oname)</div><div> oout = addfile(mask_oname,"c")</div><div> oout->odat_mask = odat_mask</div><div> else</div><div> print("Reading mask off file.")</div><div><br></div><div>;---Read the new mask from the NetCDF file</div><div> omask = addfile(mask_oname,"r")</div><div> odat_mask = omask->odat_mask</div><div> end if</div><div><br></div><div>;; ndat</div><div>if(CREATE_MASK) then</div><div> print("Creating the mask file...")</div><div><br></div><div>;---Create a new mask using a shapefile of USA</div><div> ndims = dimsizes(ndat)</div><div> optn = True</div><div> optn@return_mask = True</div><div>; optn@debug = True</div><div> optn@minlon = minlon ; Makes the shapefile masking</div><div> optn@maxlon = maxlon ; go faster.</div><div> optn@minlat = minlat</div><div> optn@maxlat = maxlat</div><div> ndat_mask = shapefile_mask_data(ndat,shp_fname,optn)</div><div><br></div><div>;---Write new mask to file</div><div> system("rm -f " + mask_nname)</div><div> nout = addfile(mask_nname,"c")</div><div> nout->ndat_mask = ndat_mask</div><div> else</div><div> print("Reading mask off file.")</div><div><br></div><div>;---Read the new mask from the NetCDF file</div><div> nmask = addfile(mask_nname,"r")</div><div> ndat_mask = nmask->ndat_mask</div><div> end if</div><div><br></div><div>;---Create masked data array</div><div> omask = where(odat_mask.eq.1,odat,odat@_FillValue)</div><div> copy_VarMeta(odat,omask)</div><div><br></div><div>;---Create masked data array</div><div> nmask = where(ndat_mask.eq.1,ndat,ndat@_FillValue)</div><div> copy_VarMeta(ndat,nmask)</div><div><br></div><div><br></div><div> ;ind32=ind(ndtooned(nmask).eq.32)</div><div> ;print(ind32)</div><div><br></div><div> dimsnmk=dimsizes(nmask)</div><div> do i=0,dimsnmk(0)-1</div><div> do j=0,dimsnmk(1)-1</div><div> if(any(ndat(i,j).eq.32)) then</div><div> print(i+" "+j+" "+" "+ndat(i,j))</div><div> end if</div><div> end do</div><div> end do</div><div><br></div><div>;;;;;;;;;;;;;;;;;;;;;;;;;;;;</div><div>;; create plot</div><div>;;;;;;;;;;;;;;;;;;;;;;;;;;;;</div><div>plot=new(2,"graphic")</div><div>info=(/ \</div><div> "1 urban_and_built-in" ,\</div><div> "2 dryland_cropland" ,\</div><div> "3 irrigated_cropland" ,\</div><div> "4 mixed_dry_cropland" ,\</div><div> "5 cropland_grassland" ,\</div><div> "6 cropland_woodland" ,\</div><div> "7 grassland" ,\</div><div> "8 shrubland" ,\</div><div> "9 mixed_shrub_grassland",\</div><div> "10 savanna" ,\</div><div> "11 deciduous_broadleaf" ,\</div><div> "12 deciduous_needleleaf" ,\</div><div> "13 evergreen_broadleaf" ,\</div><div> "14 evergreen_needleleaf" ,\</div><div> "15 mixed_forest" ,\</div><div> "16 water" ,\</div><div> "17 herbaceous_wetland" ,\</div><div> "18 wooded_wetland" ,\</div><div> "19 barren_or_sparse" ,\</div><div> "31 low_intensity_res" ,\</div><div> "32 mid_intensity_res" ,\</div><div> "33 high_intensity_res"/)</div><div><br></div><div>ninfo=dimsizes(info)</div><div><br></div><div> ninfo = dimsizes(info) ; # of classifications</div><div><br></div><div><br></div><div>;************************************************</div><div>; create plot</div><div>;************************************************</div><div> pltType = "png" ; ps, pdf, png, x11, eps</div><div><br></div><div><br></div><div><br></div><div> colorscheme =(/"red","lawngreen","mediumseagreen","springgreen",\</div><div> "limegreen","darkolivegreen3","green","lightsalmon","green3",\</div><div> "gold","palegreen","olivedrab3","chartreuse3","chartreuse", \</div><div> "darkgreen","blue","hotpink2","hotpink","light cyan", \</div><div> "khaki1","violetred1","purple3"/)</div><div><br></div><div> ncolors = dimsizes(colorscheme)</div><div> if (ninfo.ne.ncolors) then ; make sure # of colors match categories (classes)</div><div> print("size mismatch: ninfo="+ninfo+" ncolors="+ncolors)</div><div> exit</div><div> end if</div><div><br></div><div> wks = gsn_open_wks(pltType,"plot_origin_new_lc_test")</div><div> plot=new(2,graphic)</div><div> res = True ; plot mods desired</div><div> res@gsnDraw = False</div><div> res@gsnFrame = False</div><div> ;res@gsnMaximize = True ; ps, pdf</div><div> res@gsnAddCyclic = False</div><div><br></div><div> ;res@mpMinLatF = min(odat&lat)</div><div> ;res@mpMaxLatF = max(odat&lat)</div><div> ;res@mpMinLonF = min(odat&lon)</div><div> ;res@mpMaxLonF = max(odat&lon)</div><div><br></div><div> res@cnFillOn = True ; color Fill</div><div> res@cnFillMode = "RasterFill"</div><div> res@cnLinesOn = False ; Turn off contour lines</div><div> res@cnLineLabelsOn = False</div><div> res@cnLevelSelectionMode = "ExplicitLevels" ; set explict contour levels</div><div> res@cnLevels = tobyte(ispan(1,ninfo-1,1)+1)</div><div> res@cnFillPalette = colorscheme ; distinct colors for categories</div><div> res@gsnSpreadColors = False ; use each color sequentially</div><div> res@cnInfoLabelOn = False</div><div><br></div><div><br></div><div> res@mpFillOn = False</div><div><br></div><div><br></div><div> res@gsnAddCyclic = False ; regional data</div><div> ;res@trGridType = "TriangularMesh"</div><div> res@tmXBLabelFontHeightF = 0.01 ; Make lon & lat text smaller</div><div> res@tmYLLabelFontHeightF = res@tmXBLabelFontHeightF</div><div><br></div><div> res@gsnLeftString = ""</div><div> res@gsnRightString = ""</div><div><br></div><div> res1 = True</div><div> res1 = res</div><div> res1@mpMinLatF = minlat</div><div> res1@mpMaxLatF = maxlat</div><div> res1@mpMinLonF = minlon</div><div> res1@mpMaxLonF = maxlon</div><div><br></div><div> res1@tmXBMode = "Explicit"</div><div> res1@tmXBValues = (/116.0,116.1,116.2,116.3,116.4,116.5,116.6,116.7,116.8/)</div><div> 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",\</div><div> "116.6~S~o~N~E","116.7~S~o~N~E","116.8~S~o~N~E"/)</div><div><br></div><div> res1@tmYLMode = "Explicit"</div><div> res1@tmYLValues = (/39.6,39.7,39.8,39.9,40.0,40.1,40.2/)</div><div> 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"/)</div><div><br></div><div> res1@tiMainString = "confused LC"</div><div> res1@lbLabelBarOn = False</div><div> plot(0) = gsn_csm_contour_map_ce(wks, omask, res1) ; create plot</div><div><br></div><div> res2 = True</div><div> res2=res</div><div> res2@mpMinLatF = minlat</div><div> res2@mpMaxLatF = maxlat</div><div> res2@mpMinLonF = minlon</div><div> res2@mpMaxLonF = maxlon</div><div> res2@lbLabelBarOn = False</div><div><br></div><div> res2@tmXBMode = "Explicit"</div><div> res2@tmXBValues = (/116.0,116.1,116.2,116.3,116.4,116.5,116.6,116.7,116.8/)</div><div> 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",\</div><div> "116.6~S~o~N~E","116.7~S~o~N~E","116.8~S~o~N~E"/)</div><div><br></div><div> res2@tmYLMode = "Explicit"</div><div> res2@tmYLValues = (/39.6,39.7,39.8,39.9,40.0,40.1,40.2/)</div><div> 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"/)</div><div><br></div><div> res2@tiMainString = "reclassified LC"</div><div> plot(1) = gsn_csm_contour_map_ce(wks,nmask,res2)</div><div><br></div><div> pres = True</div><div> pres@gsnFrame = False</div><div> pres@gsnPanelLabelBar = True</div><div> pres@lbLabelFontHeightF = 0.01</div><div> pres@gsnPanelBottom = 0.05</div><div> pres@lbLabelPosition = "Center" ; label position</div><div> pres@lbLabelAlignment = "BoxCenters" ; label orientation</div><div> pres@cnLevels = (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,\</div><div> 31,32,33/)</div><div> pres@lbLabelStrings = (/"1","2","3","4","5","6","7","8","9","10",\</div><div> "11","12","13","14","15","16","17","18","19",\</div><div> "31","32","33"/)</div><div> ;pres@lbLabelStrings = ispan(0,ninfo,1)</div><div> ;pres@lbLabelStride = 1</div><div> pres@lbLabelAutoStride = True</div><div><br></div><div> gsn_panel(wks,plot,(/1,2/),pres)</div><div><br></div><div> rtxt = True</div><div> rtxt@txJust = "TopLeft"</div><div> rtxt@txFontHeightF = 0.0125</div><div><br></div><div> ; Add text: rows x columns of text (arbitrary)</div><div>; Usually must play with xx, yy and txFontHeightF</div><div><br></div><div> nrow = 5 ; # rows</div><div> ncol = 6 ; # columns</div><div><br></div><div> n = -1</div><div> xx = 0.015 ; arbitrary</div><div> do nc=0,ncol-1</div><div> yy = 0.20</div><div> do nr=0,nrow-1</div><div> n = n+1</div><div> if (n.le.(ninfo-1)) then</div><div> gsn_text_ndc (wks,info(n),xx,yy,rtxt)</div><div> yy = yy - 2*rtxt@txFontHeightF</div><div> end if</div><div> end do</div><div> xx = xx + 0.20</div><div> end do</div><div><br></div><div> ;draw(plot)</div><div> frame(wks)</div><div><br></div><div>end</div></div><div><br></div><div><br></div><div><br></div><div>++++++++++++++++++++++++++++++++++</div><div><br></div><div>i don't know if the nested pic can show, so i attach it </div><div><br></div><div>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.</div><div>the 31-32 exists in data.</div><div><br></div><div>can the ncl script match the landcover index automatically?</div><div><br></div><div><br></div><div>would you please give me some advices for this problem which i have been confused some days？</div><div><br></div><div>thanks</div></div>