[R-sig-Geo] How to crop a shapefile ?

[Dan Putler]

Hi,

I have used the sp pacakge, along with several other tools, to do
exactly this (although using road network files from other sources). It
can be done entirely within R via sp, and I've written R only code to do
it, but found it very slow (someone else on the list may have an idea as
to how to make this more efficient).

Let me describe the method I use. One question I have is do you want the
road segments to fit within a bounding box, or do you want them to fit
into some other administrative area? If you have a polygon file of the
administrative area, you can use it to crop things even more precisely.

The steps to the approach follow:

1. Extract the the first ("from") and last ("to") point of each polyline
and append these two points to the SpatialLineDataFrame's data frame.
(Note: This is what I find to be slow doing in R via the sp package.
What I use to do this quickly is a C program I wrote using the shapelib
library. My guess is that one of the sp package authors can suggest ways
to speed this up).

2. Test to see whether either of the from or to points of each polyline
are within the bounding box of interest.

3. Create a new SpatialLineDataFrame object based on selecting only
those polylines that meet the bounding box test.

This approach will likely have some road segments that extend beyond the
bounding box, but have one of its endpoints in the bounding box. Most
road network files snap all polylines at intersections, but this is not
true with OSM "ways" (some of which are extremely long). As a result,
some of the "fuzz" polylines may extend well outside the bounding box.
In addition, long ways that run through the bounding box, but neither
begin or end within it, will not be selected.

Dan

On Thu, 2009-12-17 at 11:07 +0100, Baptiste Coulmont wrote:
> Hi,
> 
> I'm working with R and various shapefiles created by the OpenStreetMap
> project. Those shapefiles (the roads in particular) are quite big.
> 
> How can I crop (or "trim", or "prune" those objects -- once in a
> SpatialLineDataFrame -- if I need to work not on a national level but
> on a municipal level ?
> 
> Is it possible to create a subset of a smaller zone ?
> Coordinates:
>          min       max
> r1  2.560776  6.403847
> r2 49.493613 51.505419
> 
> Thansk,
> B. Coulmont
> 
-- 
Dan Putler
Sauder School of Business
University of British Columbia