Re: Subset of a shapefile

Thank you Mary, and thank you ncl-talk. What service!

Jim

On 4/6/2013 22:46, Mary Haley wrote:
> Hi Jim,
>
> You need to copy the gsn_add_shapefile_polylines function
> from $NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl and make some
> minor modifications to it, since it wants to draw everything in the
> shapefile.
>
> I've created an example, since this is likely a highly requested
> feature. See example shapefiles_12.ncl at:
>
> http://www.ncl.ucar.edu/Applications/shapefiles.shtml#ex12
>
> This code:
>
> ;---Attach specified highways as polylines on the map.
> highways = (/"I- 5","I- 82","I- 90"/)
> pres = True
> pres@gsLineColor = "blue"
> pres@gsLineThicknessF = 3.0 ; default is 1.0
> poly = gsn_add_shapefile_polylines_subset(wks,map,shp_filename,"ROUTE",highways,pres)
>
> is where you specify which features from the shapefile you want to
> draw, and what the name of the variable is on the shapefile that
> contains these features ("ROUTE").
>
> Note to others: this script easily be used for any polyline shapefile
> that you want to draw only a subset of features from. You can also
> modify this code to indicate a list of features that you *don't* want
> to draw. Simply remove the operator ".not" from both instances of
> this "if" statement:
>
> if(.not.any(feature_names(i).eq.vlist)) then
>
>
> --Mary
>
>
> On Apr 5, 2013, at 6:32 PM, Jim Means wrote:
>
>> Hello, I am trying to plot a single interstate highway (I-5) in NCL.
>> I have downloaded a shapefile of interstates and other roads from
>> NOAA
>> (http://www.nws.noaa.gov/geodata/catalog/transportation/html/interst.htm),
>> and by stealing code from one of the NCL examples I'm able to load
>> the shapefile into NCL and plot the road. Unfortunately there are
>> other roads that I don't want to plot that show up also. So far, I've
>> been unable to figure out how to just plot a single road. My code
>> looks like this:
>>
>> f = addfile("in101503.shp", "r") ; Open shapefile
>>
>> ;
>> ; Read data off shapefile
>> ;
>> segments = f->segments
>> geometry = f->geometry
>> segsDims = dimsizes(segments)
>> geomDims = dimsizes(geometry)
>>
>> numFeatures = geomDims(0)
>> ;
>> ; 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)
>>
>> interstates = str_match_ind_ic(f->ROUTE, "I- 5")
>>
>> The last statement gives me the indices that correspond to I-5 in the
>> shapefile. I think I should be able to use this to just plot that
>> single road, but after much experimentation and trying to follow the
>> examples, I haven't been able to do it. Any hints would be greatly
>> appreciated.
>>
>> Thank you,
>>
>> Jim Means
>>
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
Received on Sat Apr 6 16:52:20 2013