5

I'm attempting to perform a gIntersection of a big grid on a shapefile for country boundaries, so that I can calculate the % area of each grid cell that lies within the country. It takes a few seconds for smaller countries, but takes ages for larger ones.

I have a method that cut down the time substantially, but it still takes an unreasonable amount of time for larger countries. Here's a reproducible example, that uses gIntersects to subset the grid, which cut down on time substantially.

#### demo script for stackoverflow ####
# set up demo grid
lat <- seq(from = -11, to = 55, by = 0.5)
long <- seq(from = 70, to = 150, by = 2/3)
demo.grid <- SpatialPoints(coords = cbind(rep(long, times = length(lat)), rep(lat, each = length(long))))
demo.grid <- SpatialPixels(demo.grid, tolerance = 0.000001, proj4string = CRS("+proj=longlat +datum=WGS84"))
demo.grid.poly <- as.SpatialPolygons.GridTopology(getGridTopology(demo.grid), proj4string = CRS("+proj=longlat +datum=WGS84"))

# looks like a grid
plot(demo.grid.poly)

# get my file - I used the SpatialPolygonsDataFrame for level 0 administrative boundaries for China, downloaded from www.gadm.org
CHN <- get(load(file.choose()))
plot(CHN, add = TRUE, border = "blue") # China is located in the right place

CHN.p <- spTransform(CHN, CRS("+proj=aea +lat_1=7 +lat_2=-32 +lat_0=-15 +lon_0=125 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))
demo.grid.poly.p <- spTransform(demo.grid.poly, CRS("+proj=aea +lat_1=7 +lat_2=-32 +lat_0=-15 +lon_0=125 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))

plot(demo.grid.poly.p)
plot(CHN.p, add = TRUE, border = "blue")

# looks good

system.time(subset <- gIntersects(CHN.p, demo.grid.poly.p, byid = TRUE)) # find the grid cells that intersect with China
summary(subset)
demo.grid.poly.subset <- demo.grid.poly.p[subset] # make new set of gridcells that are subsetted to only those that intersect with China

plot(demo.grid.poly.subset)
plot(CHN.p, add = TRUE, border = "blue")
# still looks good
system.time(intersection <- gIntersection(CHN, demo.grid.poly.subset, byid = TRUE)) # create intersection SPDF. Go get sandwich, sleep, exercise, take a small vacation. Learn piano.

This is related to a question I asked a while ago (https://stackoverflow.com/questions/16918767/gintersection-on-very-large-spatial-objects). I'm facing a similar issue with a new dataset, and have a new method for solving the problem that has a new problem.

  • 1
    Isn't there a gSimplify method or similar that reduces the fine detail in the edges? Should speed things up considerably. – Ari B. Friedman Nov 26 '13 at 23:55
  • Sorry it took me so long to respond! It looks like gSimplify will be the winner. I got the process started, we will see how long it takes. It got through 4% in about 10 minutes, so I bet this will work if I leave it overnight. Thanks! – LightonGlass Jan 8 '14 at 18:52
  • Thanks for coming back to confirm what worked. I posted my comment as an answer since it seemed to work. – Ari B. Friedman Jan 8 '14 at 19:02
1

You could try a two step process:

  1. Use gSimplify to reduce the complexity of your polygons' edges.
  2. Then run gIntersection as normal.
2

An alternative would be to use the raster package. You can rasterize a polygon, mask to the polygon raster and then calculate information on the subset raster. This is quite fast compared to what you currently have implemented.

require(raster)
require(sp)

# create example data
r <- raster(ncol=500, nrow=500)
  r[] <- runif(ncell(r),0,1)
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), 
                              Polygons(list(Polygon(cds2)), 2)))

plot(r)
  plot(polys, add=TRUE)         

# Get first polygon 
f <- polys[1,]

# Crop raster to polygon. This is used as reference raster in rasterize 
cr <- crop(r, extent(f), snap="out")                    

# Rasterize polygon  
fr <- rasterize(f, cr)

# Mask to polygon
r.sub <- mask(x=cr, mask=fr)

plot(r.sub, main="Subset raster")
  plot(f, add=TRUE) 
  • This looks like it may work in some cases, but not in mine. It looks like rasterize() takes values from the center of a grid cell, but I need to get the % area. If I had a grid cell the center of which was in, say, Pennsylvania, but very close to the border with Ohio, so it was 40% in Ohio and 60% in PA, it looks like rasterize() would say that cell was 100% PA. I need the 40/60 split. – LightonGlass Dec 16 '13 at 19:30
  • Then you may need to get out of the raster topology class entirely and treat your raster cells as true polygon topology in the gIntersects function which is meant to deal with vectors. You can coerce to a sp polygon object using: as(x, "spatialPolygonDataFrame"). Using a polygon object may speed things up but I would not hold my breath. I have also noticed that in the current release of sp the over function have seen notable speed gains. This size problem may be difficult in R. – Jeffrey Evans Dec 16 '13 at 20:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.