Subset of a shapefile

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 Fri Apr 5 18:32:47 2013