1

I have two huge shapefiles. One contains polygons of one kind of vegetation and the other one representing a grid. A situation as follows but, of course, much bigger:

enter image description here

I need:

  1. the efficient way to open and handle both files using R, because my computer "suffers" opening both using readOGR or shapefile even if I have a 32 GB RAM i7;
  2. I need to get a matrix-like (or dataframe) containing the number of polygons for each square within the grid, something like this (from the previous example):

enter image description here

Improve this question
6
  • 2
    How big is a big shapefile? Is it big because it has a lot of features, or because it has a lot of attributes, or because it has a lot of detail in its features? How big in mega/kilo/giga bytes are the various shapefile parts? Can you read them in at all? Have you tried st_read from the sf package? Is the grid completely regular and axis-aligned? If it is then you don't need to read it in, you can work out things from the geometry only (extent and resolution).
    – Spacedman
    Commented Feb 18, 2019 at 23:44
  • The grid is not that big: just 9700 KB of perfectly 4km * 4 km pixels (a little bit more than 73000 pixels). The other one has more than 550000 patches and weighs almost 1300000 KB. Yes, I tried st_read but failed using it to get the reckon required.
    – perep1972
    Commented Feb 19, 2019 at 0:04
  • 1
    I'd try expressing your grid as a raster - i.e. r <- raster(extent(xmin, xmax, ymin, ymax), nrows = nr, ncols = nc) and then raster::rasterize(polygons, r, fun = 'count') - but you'll need an up to date version of raster, and definitely experiment with a subset of the patches to get a sense of the time required - you say you can't read the shapefile with sf::read_sf? Definitely need a solution for that first.
    – mdsumner
    Commented Feb 19, 2019 at 0:24
  • 1
    Ah, what I said can't work - the lookup is very approximate, just using the centre of each cell.
    – mdsumner
    Commented Feb 19, 2019 at 0:51
  • 1
    If you can't read your shapefile entirely then have you tried using an SQL clause to read it in chunks? Or to split it into smaller shapefiles? Or convert to GeoPackage and read in chunks?
    – Spacedman
    Commented Feb 19, 2019 at 8:10

1 Answer 1

Reset to default
4

You can intersect your grid with the polygons and aggregate on the number of objects per grid square.

read_sf is a lot faster than readOGR but with a much larger dataset you may still struggle with the spatial operations - if this is a problem then you could try PostGIS.

# Read in all the regions in South America (282 polygons)
sam <- read_sf(dsn = "~/workspace/_temp/gridsmall.gpkg", layer = "samerica")

The northern tip:

# Read the 1° grid covering all of S.Am (4278 grid squares) 
samgrid <- read_sf(dsn = "~/workspace/_temp/gridsmall.gpkg", layer = "samgrid")


# Intersect the two, which cuts all the polygons by each grid boundary:
sam.int <- st_intersection(sam, samgrid)

Each grid had a unique id number, 'id' - use this to summarise the number of names (from the region name of each polygon) in each grid square:

sam.agg <- aggregate(id ~ name, data = sam.int, FUN = length)

This will give a data frame with the number of polygons per grid square.

Improve this answer
2
  • Thanks a lot, @Simbamangu. Was a little bit more complicated than that, but you gave me a great tip. Thank you. Best!
    – perep1972
    Commented Feb 19, 2019 at 12:59
  • 1
    @perep1972, you're welcome - suggest you look into PostGIS as well (I did this operation on admin level 1 regions worldwide (5k polygons), intersecting with 1 degree grid squares (54k cells) and it took ~ 5 mins). Also, please click the check mark next to the answer to mark it 'answered'.
    – Simbamangu
    Commented Feb 19, 2019 at 13:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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