Piping OGR/GDAL in R

Question

I have to read in a small portion of a huge ESRI shapefile in R. I am doing this in two steps:

Step 1: I use ogr2ogr to clip the shapefile to my bounding box:

ogr2ogr -clipsrc xMin yMin xMax yMax outfile.shp infile.shp

Step 2: I read this into R with rgdal:

df = readOGR(dsn="/path", layer="outfile")

The problem is that I have to do this for multiple files, and it's hard to keep track of the OGR operations which generated each of these separate files. Is there any way to pipe ogr2ogr in R, so that step 1 is done on the fly?

[Here's an example](http://stackoverflow.com/questions/13982773/crop-for-spatialpolygonsdataframe/13986029#13986029) of how to clip a polygon using `rgeos::gIntersection`. (It's basically the same as @mdsumner's answer below, but with a reproducible example and a purty picture.) — Josh O'Brien, Sep 10 '13 at 15:04

score 5 · Accepted Answer · answered Sep 10 '13 at 09:40

5

Try using the system call. Hard to tell without code and data, but you should be able to create a list of either shapefiles to process (if you have multiple shapefiles) or coordinates to process if you have multiple bounding boxes you want to clip from the shapefile.

work.dir <- "C:/Program Files (x86)/FWTools2.4.7/bin" # use your FWTools location
setwd(work.dir)
out.shape.file <- "foo2.shp"
in.shape.file <- "foo1.shp"
x.min <- 100
y.min <- 100
x.max <- 200
y.max <- 200
system(paste("ogr2ogr -clipsrc", x.min, y.min, x.max, y.max, 
      out.shape.file, in.shape.file))

answered Sep 10 '13 at 09:40

SlowLearner

7,907
11
49
80

Worked like a charm. Many thanks, SlowLearner. (PS: If this is slow, I wonder what you'll be like when you pick up speed.) – user702432 Sep 10 '13 at 10:00
Well, it's all relative isn't it - we both have a long way to go compared to the R veterans of Stack Overflow... – SlowLearner Sep 10 '13 at 10:02
True for me, at least. – user702432 Sep 10 '13 at 10:08

score 1 · Answer 2 · answered Sep 10 '13 at 10:39

Note that you can clip with the rgeos package, something like this (I use raster to make creating a simple mask polygon easy):

library(rgeos)
library(raster)
library(rgdal)
## as per your example
df = readOGR(dsn="/path", layer="outfile")
## create a simple polygon layer and intersect
res <- gIntersection(df, as(extent(x.min, x.max, ymin, y.max), "SpatialPolygons"), byid = TRUE)

That might not work exactly as is, but could be worth exploring if reading the entire shapefile is not a problem.

Hi mdsumner... The problem is that the files are too large to load directly in R. That's why I was doing it in two steps. SlowLearner's solution works beautifully. — user702432, Sep 11 '13 at 06:48

Piping OGR/GDAL in R

2 Answers2