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.