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

Mary Haley haley at ucar.edu
Mon Oct 26 09:48:38 MDT 2015


I think these are map outlines that you see.  Try setting:

res at mpOutlineOn = False

--Mary


On Mon, Oct 26, 2015 at 9:21 AM, dyjbean at 163.com <dyjbean at 163.com> wrote:

> 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
>
> _______________________________________________
> 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/20151026/7e4d1cee/attachment-0001.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: bg.jpg
Type: image/jpeg
Size: 20365 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151026/7e4d1cee/attachment-0001.jpg 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: InsertPic_.png
Type: image/png
Size: 247761 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151026/7e4d1cee/attachment-0001.png 


More information about the ncl-talk mailing list