5

I have a SpatialPolygonsDataFrame that I created by reading in a shapefile using readOGR in the rgdal package. I am attempting to use it to generate a sampling grid using spsample in the sp package to be used for interpolation from survey data collected in the area. However, the SpatialPointsDataFrame encompasses a much larger area than the survey, and as a result, the interpolation is predicting values far away from where any surveys were conducted. The survey data and the shapefile were both projected using the same proj4string.

I would like to subset the SpatialPolygonsDataFrame using the coordinates set by the survey stations, but I'm not sure where the relevant values are stored in the object.

I'm afraid that I can't provide the relevant data as the shapefile is not hosted online. However, I will borrow some code from Paul Hiemstra's response to this post for the Netherlands:

library(ggplot2)
library(sp)
library(automap)
library(rgdal)

#get the spatial data for the Netherlands
con <- url("http://gadm.org/data/rda/NLD_adm0.RData")
print(load(con))
close(con)

class(gadm)
bbox(gadm)
>         min      max
>r1  3.360782  7.29271
>r2 50.755165 53.55458

Let's say that the surveys were conducted in this area:

bbox(surveys)
    >         min      max
    >r1     4.000    7.000
    >r2    51.000   53.000

How can I crop out that area of the SpatialPolygonsDataFrame?

EDIT: This question appears to answer mine. Apologies for not searching hard enough (although the comments did give me some sense of where to turn with rgeos). However, gIntersection causes R to crash...

Community
  • 1
  • 1
jslefche
  • 4,379
  • 7
  • 39
  • 50
  • The short answer is probably 'package:rgeos'. Sorry I dont have time for the long answer right now. Hence its only a comment... – Spacedman Apr 18 '12 at 16:42
  • 1
    Look at the gIntersection() function in the package rgeos – Sophia Apr 18 '12 at 16:49
  • My vague recollection of this is that one of the spatial packages has an overlay function that returns a boolean of the polygons that match. You can then use taRifx::subsetSPDF to grab only those polygons. – Ari B. Friedman Apr 18 '12 at 18:24
  • 1
    It's unclear from the question whether you want to crop the polygons to the bounding box of the points (there's only one object in your example gadm), or whether you want to remove any objects that do not intersect the points. The first is easy using over() and "[" in sp, the second requires topology tools and rgeos is the go, but do you really want to cut out the bounding box of the points? – mdsumner Apr 18 '12 at 20:49
  • Sorry, should have said "the second is easy using over()", "the first requires topology" . . . – mdsumner Apr 18 '12 at 22:56
  • @Sophia thanks for the lead. I am trying `gIntersection` currently, but it is taking quite a while. Not sure if that is par for the course when working with shapefiles... – jslefche Apr 19 '12 at 16:45

1 Answers1

-1

Depending upon the size of the polygons you could do something like

range = cbind(c(4,7), c(51,53))
centroids <- coordinates(spdf)

spdf.subset <- spdf[centroids[,1] > range[1,1] &
                    centroids[,1] < range[2,1] &
                    centroids[,2] > range[1,2] &
                    centroids[,2] < range[2,2],]
fgregg
  • 3,173
  • 30
  • 37