Re: Subset of a shapefile

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

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sat Apr 6 16:46:39 2013