[ncl-talk] extra profile in plot with shapefile, how can i get rid of them?

dyjbean at 163.com dyjbean at 163.com
Mon Oct 26 09:21:11 MDT 2015


hi,
  i made a plot with shapefile , but there exist some extra lines in figure, as follows:
    

 there some extra lines on figure above. how can i get rid of them?


the following is my 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"
forgname="beijing_globecover_liuhuan30_usgs_changed.txt"
;fconname="../../globecover_huabei_liuhuan_usgs_new_allocate1.txt"
fconname="beijing_globecover_liuhuan30_usgs_changed_lir-mir-hir.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"

;; 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"

;
; 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 = "jingjinji.shp"

;---Name of file to write mask to or to read mask from.
   mask_oname = "jingjinji_mask_300m.nc"
   mask_nname = "jingjinji_mask_1km.nc"

;---Rough area we are interested in. Everything outside this will be masked.
   minlon = 113.
   maxlon = 120.
   minlat =  36.0
   maxlat =  43.0

;; 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)


print(odat(onlat-1,omlon-1))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; 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","SpringGreen2","SpringGreen3","springgreen",\
       "limegreen","darkolivegreen3","green","lightsalmon","green3",\
        "LightSeaGreen","palegreen","olivedrab3","chartreuse2","chartreuse3",       \
        "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_lc_jingjinji_new")
  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((/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)
  res at cnFillPalette    = colorscheme        ; distinct colors for categories
  res at gsnSpreadColors  = False              ; use each color sequentially
  res at cnInfoLabelOn   = False

  ;res at lbLabelPosition  = "Center"           ; label position
  ;res at lbLabelAlignment = "BoxCenters"       ; label orientation
  ;res at lbLabelStrings    = (/"1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16",\
                        ;"17","18","19","31","32","33"/)
 ;res at lbLabelStrings   = ispan(0,ninfo,1)
  ;res at lbLabelStride    = 1
  ;res at lbLabelAutoStride = False            ; in V6.1.0 and up, this is defaulted
                                            ; to True, and overrides lbLabelStride
                                            ; default. Set this to ensure that
                                            ; your label scheme is preserved.
  ;res at pmLabelBarHeightF        =  0.075
  ;res at pmLabelBarWidthF         =  0.60      ; default is 0.6
  ;res at pmLabelBarOrthogonalPosF =  0.08      ; move up smidge

  res at mpFillOn         = False

  ;res at mpGridAndLimbOn   = True
  ;res at mpGridLatSpacingF = 5.0
  ;res at mpGridLonSpacingF = 5.0
  ;res at mpGridLineDashPattern = "2"

  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   = (/113.0,114.0,115.0,116.0,117.0,118.0,119.0,120.0/)
  res1 at tmXBLabels  =  (/"113.0~S~o~N~E","114.0~S~o~N~E","115.0~S~o~N~E","116.0~S~o~N~E","117.0~S~o~N~E",\
                        "118.0~S~o~N~E","119.0~S~o~N~E","120.~S~o~N~E"/)


  res1 at tmYLMode    = "Explicit"
  res1 at tmYLValues   = (/36.0,37.0,38.0,39.0,40.0,41.0,42.0,43.0/)
  res1 at tmYLLabels  = (/"36.0~S~o~N~N","37.0~S~o~N~N","38.0~S~o~N~N","39.0~S~o~N~N","40.0~S~o~N~N",\
                       "41.0~S~o~N~N","42.0~S~o~N~N","43.0~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 tmXBMode    = "Explicit"
  res2 at tmXBValues   = (/113.0,114.0,115.0,116.0,117.0,118.0,119.0,120.0/)
  res2 at tmXBLabels  =  (/"113.0~S~o~N~E","114.0~S~o~N~E","115.0~S~o~N~E","116.0~S~o~N~E","117.0~S~o~N~E",\
                        "118.0~S~o~N~E","119.0~S~o~N~E","120.~S~o~N~E"/)

  res2 at tmYLMode    = "Explicit"
  res2 at tmYLValues   = (/36.0,37.0,38.0,39.0,40.0,41.0,42.0,43.0/)
  res2 at tmYLLabels  = (/"36.0~S~o~N~N","37.0~S~o~N~N","38.0~S~o~N~N","39.0~S~o~N~N","40.0~S~o~N~N",\
                       "41.0~S~o~N~N","42.0~S~o~N~N","43.0~S~o~N~N"/)

  res2 at lbLabelBarOn   =  False
  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.015
  pres at lbLabelPosition  = "Center"           ; label position
  pres at lbLabelAlignment = "BoxCenters"       ; label orientation
  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 = False

  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
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

what kind of properties should i add in the script? 
or something is wrong with the shapefile?


thanks for your advices





dyjbean at 163.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151026/1d32335b/attachment-0001.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/png
Size: 247761 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151026/1d32335b/attachment-0001.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/jpeg
Size: 20365 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151026/1d32335b/attachment-0001.jpe 


More information about the ncl-talk mailing list