[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