[ncl-talk] Most data disappeared after mask data with shapefile according to the example on the website.
Rick Brownrigg
brownrig at ucar.edu
Wed Aug 10 14:02:51 MDT 2016
Hi,
I looked at this, but can not really tell what's going on without knowing
more about the data. There are two shapefiles, shaanxi_city_l.shp, and
xian.shp. I presume one of these is used to draw the subregional
boundaries in your plot (which one)? What is the other? Are these both
polygonal shapefiles? Perhaps plot them in two different colors to make
sure they are what you think they are?
Without access to the data, I'm not sure what else to suggest....
Wish I could be of more help.
Rick
On Wed, Aug 10, 2016 at 2:54 AM, grace <313695096 at qq.com> wrote:
> Hi:
> All,I am trying to mask data with shapefile according to the example
> shapefiles_14_mask.ncl on the NCL website.
> But most data disappeared after mask.
> The right plot is the plot after mask,the left one is the original plot.
>
> This is my script:
> ; Example script to produce precip plots with satationdata,
>
> 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"/public/home/huanglei/map/shapefile_mask_data.ncl"
> begin
> station_file ="/public/home/huanglei/data/20160724/stationdata_20160724.
> txt"
> data = asciiread(station_file,-1,"string")
> num_att = numAsciiCol(station_file) ;return the number of columns of the
> stationfile
> print(num_att)
>
> lon_station1 = stringtofloat(str_get_field(data(0:1889),2," "))
> lat_station1 = stringtofloat(str_get_field(data(0:1889),3," "))
> obs_precip1 = stringtofloat(str_get_field(data(0:1889),4," "))
>
>
> olon = fspan(105.5,110.9,100)
> olat = fspan(31.5,35.8,100)
>
> olon!0 = "lon"
> olon at long_name = "lon"
> olon at units = "degrees-east"
> olon&lon = olon
> olat!0 = "lat"
> olat at long_name = "lat"
> olat at units = "degrees_north"
> olat&lat = olat
> print("......")
> obs_precip1 at _FillValue = -999.000000
> rscan = (/0.5,0.2,0.1/)
> griddata = obj_anal_ic_deprecated(lon_station1,lat_station1,obs_
> precip1,olon,olat,rscan,False)
> ;;;;;;;;;;;;;;;;;;;;;;mask data with shapfile;;;;;;;;;;;;;;;;;;;
> shp_filename = "/public/home/huanglei/map/shaanxi_city_l.shp"
> final= shapefile_mask_data(griddata,shp_filename,True)
> ; exit
> wks = gsn_open_wks ("pdf","rain20160724")
>
> res = True
> res at gsnMaximize = True
> res at gsnDraw = False
> res at gsnFrame = False
>
> ;>--------------------------------------------<
> ; set for the map
> ;>--------------------------------------------<
> res at mpMinLatF = 31.5
> res at mpMaxLatF = 35.8
> res at mpMinLonF = 105.5
> res at mpMaxLonF = 110.9
> res at tmXBMode = "Explicit"
> res at tmXBValues = (/106,107,108,109,110/)
> res at tmXBLabels = (/"106~S~o~N~E","107~S~o~N~E",
> "108~S~o~N~E","109~S~o~N~E","110~S~o~N~E"/)
> ; res at tmXBMinorValues = fspan(97,109,30)
> ; res at tmXBMinorOn = True
>
> res at tmYLMode = "Explicit"
> res at tmYLValues = (/32,33,34,35/)
> res at tmYLLabels = (/"32~S~o~N~N","33~S~o~N~N","34~S~o~N~N","35~S~o~N~N"/)
>
> ;res at tmYLMinorValues = fspan(26,36,20)
> ;res at tmYLMinorOn = True
>
> res at mpFillOn = True
> res at mpOutlineOn = False
> res at cnFillDrawOrder = "PreDraw"
> res at mpDataBaseVersion = "MediumRes"
> res at mpDataSetName = "Earth..4"
> res at mpAreaMaskingOn = True
> res at mpLandFillColor = "white"
> res at mpInlandWaterFillColor = "white"
> res at mpOceanFillColor = "white"
> res at mpOutlineBoundarySets = "NoBoundaries"
>
> ;>--------------------------------------------<
> ; set for the plot
> ;>--------------------------------------------<
>
> res at cnFillOn = True
> res at cnLinesOn = False
> ; res at cnLevelSpacingF = 0.2
> res at gsnSpreadColors = True
> res at lbLabelAutoStride = True
> res at gsnAddCyclic = False
>
> res at cnLevelSelectionMode = "ExplicitLevels"
> res at cnLevels = (/ .1, .4, .8, 1.6, 3.2, 6.4, \
> 12.8, 25.6, 51.2,100/)
> res at cnFillColors = (/"White","DarkOliveGreen1", \
> "DarkOliveGreen3","Chartreuse", \
> "Chartreuse3","Green","ForestGreen",
> \
> "Yellow","Orange","Red","Violet"/)
> map = gsn_csm_contour_map(wks,final,res)
>
> ;>------------------------------------------------------------<
> ; add map
> ;>------------------------------------------------------------<
>
> shp_name1 = "/public/home/huanglei/map/xian.shp"
> shp_name2 = "/public/home/huanglei/map/shaanxi_city_l.shp"
> lnres = True
> lnres at gsLineColor = "gray25"
> lnres at gsLineThicknessF = 1.5
>
> id = gsn_add_shapefile_polylines(wks,map,shp_name1,lnres)
> id2= gsn_add_shapefile_polylines(wks,map,shp_name2,lnres)
> ;>------------------------------------------------------------<
>
> ; station = asciiread("/public/home/
> huanglei/data/prcdata.txt",(/160,3/),"float")
>
> res2 = True
> res2 at gsMarkerIndex = 16
> res2 at gsMarkerSizeF = 6.
> res2 at gsMarkerColor = "Violet"
> res2 at tfPolyDrawOrder = "PostDraw"
> res2 at cnFillDrawOrder = "PostDraw"
> ;
> plots1=gsn_add_polymarker(wks,map,108.9373,34.2158,res2)
> plots2=gsn_add_polymarker(wks,map,108.9133,34.2303,res2)
> plots3=gsn_add_polymarker(wks,map,108.8908,34.1953,res2)
>
> ; plots=gsn_add_polymarker(wks,map,station(:,1),station(:,2),res2)
> delete(res2)
> draw(map)
> frame(wks)
> end
> How can I slove the problem?
>
> _______________________________________________
> 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/20160810/cb7b0134/attachment.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: FCC3B4E7 at ABD5F612.C3EBAA57.jpg
Type: image/jpeg
Size: 206822 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160810/cb7b0134/attachment.jpg
More information about the ncl-talk
mailing list