Re: ploygon shapefile subset

From: Mary Haley <haley_at_nyahnyahspammersnyahnyah>
Date: Sun Apr 07 2013 - 20:58:29 MDT

Hi Jennifer,

I think you have a typo in your program. Inside your loop across "numFeatures" you want to attach the polygons using gsn_add_polyon, and not gsn_add_shapefile_polygons.

Trying changing this line:

        poly(segNum) = gsn_add_shapefile_polygons(wks, map, "dom_2009_krig.shp", pres)

to:

        poly(segNum) = gsn_add_shapefile_polygons(wks, map, lon(startPT:endPT),lat(startPT:endPT),pres)

--Mary

On Apr 7, 2013, at 9:25 AM, Jennifer McCabe wrote:

> Hi all,
> I've been trying to fill ploygons in a shapefile different colors depending on a range of values in the "gr" column of the attribute table of the shapefile. Admittedly I am very new to NCL. I've tried using examples from the website and the new one Mary just posted yesterday for subsetting polylines. I keep getting the bellow error.
> Dimension sizes on right hand side of assignment do not match dimension sizes of left hand side
> fatal:["Execute.c":8128]:Execute: Error occurred at or near line 120 in file project2.ncl
>
> If I comment out the section causing the error I get it to draw but the entire shapefile is filled in with blue (the first color).
> Here is my code...
>
> in_grib1 = addfile("batch-mars-web240-20130319140652-TnMhMt.grib","r")
> in_grib2 = addfile("batch-mars-web240-20130319145852-7dOvP2.grib","r")
> in_grib3 = addfile("batch-mars-web240-20130319174030-8SXgmI.grib","r")
>
> function create_map(wks,BBS)
> local a, res
> begin
>
> res = True
>
> res@gsnDraw = False ; don't draw yet
> res@gsnFrame = False ; don't advance frame yet
> res@gsnMaximize = True ; maximize plot in frame
>
> res@mpProjection = "LambertConformal" ; choose projection
> res@mpLambertParallel1F = 33 ; first parallel
> res@mpLambertParallel2F = 45 ; second parallel
> res@mpLambertMeridianF = -98 ; meridian
>
> res@mpLimitMode = "Corners" ; corner method of zoom
> res@mpLeftCornerLatF = 30 ; left corner
> res@mpLeftCornerLonF = -125 ; left corner
> res@mpRightCornerLatF = 65 ; right corner
> res@mpRightCornerLonF = -35 ; right corner
>
>
> res@pmTickMarkDisplayMode = "Always" ; turn on tickmarks
>
> res@tiMainString = "BBS data"
>
> ;---Create map.
> map = gsn_csm_map(wks,res)
>
> return(map)
> end
>
> ;*************************************************
> ; main code
> ;*************************************************
> begin
>
> wks = gsn_open_wks("pdf","project2")
>
> ;---Create the map
> map = create_map(wks,"2009")
>
> ;---Attach two sets of polylines to the same map.
> lnres = True
> lnres@gsLineThicknessF = 0.5
> lnres@gsLineOpacityF = 0.5
> lnres@gsLineColor = "gray55"
> poly1 = gsn_add_shapefile_polylines(wks,map,"s_06se12_Project.shp",lnres)
>
> lnres@gsLineColor = "gray55"
> lnres@gsLineThicknessF = 0.5
> lnres@gsLineOpacityF = 0.5
> poly0 = gsn_add_shapefile_polylines(wks,map,"PROVINCE.shp",lnres)
>
> pres = True
> pres@gsFillColor = "Transparent"
> pres@gsEdgesOn = True
> pres@gsLineThicknessF = 0.06
> pres@gsLineOpacityF = 0.5
>
> f = addfile("dom_2009_krig.shp","r")
> gr = f->GRIDCODE
> lon = f->lon
> lat = f->lat
> ; Read data off shapefile
> ;
> segments = f->segments
> geometry = f->geometry
> segsDims = dimsizes(segments)
> geomDims = dimsizes(geometry)
> ;
> ; Read global attributes
> ;
> geom_segIndex = f@geom_segIndex
> geom_numSegs = f@geom_numSegs
> segs_xyzIndex = f@segs_xyzIndex
> segs_numPnts = f@segs_numPnts
> numFeatures = geomDims(0)
>
> colors = (/"blue","green","yellow","red"/)
> npoly = sum(geometry(:,geom_numSegs))
> poly = new(npoly,graphic) ; Array to hold polygons
> segNum = 0
>
> do i=0, numFeatures-1
>
> if (gr(i).ge.290 .and. gr(i).lt.400) then
> pres@gsFillColor = colors(0)
> end if
>
> if (gr(i).ge.482 .and. gr(i).lt.552) then
> pres@gsFillColor = colors(1)
> end if
>
> if (gr(i).ge.600 .and. gr(i).lt.710) then
> pres@gsFillColor = colors(2)
> end if
>
> if (gr(i).ge.720 .and. gr(i).lt.857) then
> pres@gsFillColor = colors(3)
> end if
>
> startSegment = geometry(i, geom_segIndex)
> numSegments = geometry(i, geom_numSegs)
> do seg=startSegment, startSegment+numSegments-1
> startPT = segments(seg, segs_xyzIndex)
> endPT = startPT + segments(seg, segs_numPnts) - 1
> poly(segNum) = gsn_add_shapefile_polygons(wks, map, "dom_2009_krig.shp", pres)
> segNum = segNum + 1
> end do
> end do
>
>
> poly2 = gsn_add_shapefile_polygons(wks,map,"dom_2009_krig.shp",pres)
>
>
>
> ;---Drawing the map will also draw the attached polygons.
> draw(map)
> frame(wks)
> end
>
>
> --
> Jennifer McCabe
> PhD Student, Ecology and Environmental Science
> 204 Roger Clapp Greenhouse, UMaine
> Orono, ME 04469
> p:307-413-5512
> jennifer.mccabe@maine.edu
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sun Apr 7 20:58:43 2013

This archive was generated by hypermail 2.1.8 : Sun Apr 07 2013 - 21:13:30 MDT