I want to produce a gridded map where each grid cell represents the sum of points falling within the grid using R. I have seen similar questions for R here, here, and here for Phyton. But could not get what I needed from none of them.

I have got the gridded map plotted onto a shapefile of the Neotropical region from here https://figshare.com/articles/dataset/Neotropical_region_a_shapefile_of_Morrone_s_2014_biogeographical_regionalisation/3569361 and the distribution of species on the map (each point represents an individual of each species in a specific geographic location). The distribution of points was downloaded from GBIF using the gbif R function from the dismo R package.

To get what I need, these points should be summarised by species per cell so that I get the number of species (irrespective of the number of individuals) that fall within each grid cell. Thus, I can plot the gridded map over the shapefile, coloured by the number of species.

I am a newbie to spatial analysis and have sparse knowledge about doing this in R. So my apologies if it is a very simple question. Any hint is very welcome!

Here is the code I used to get what I have so far:

# Load required packages

#retrieve occurrences of Allagoptera species from GBIF
alag <- gbif("Allagoptera", "*", geo=TRUE)

#filter those with latitude and longitude data
alagl <- subset(alag, !is.na(lon) & !is.na(lat))

#I followed the tutorial in the RFunctions blog to create the grids:
#I just replaced the shapefile with the one for the Neotropics and modified the comments accordingly. 

# Load the Neotropical shapefile.
#Data come from here:

#Function 'download.shapefile' comes from here: 


neo <- readOGR(choose.files(), "Lowenberg_Neto_2014") 

# Plot the shape to see if everything is fine.

# Create an empty raster.
grid <- raster(extent(neo))

# Choose its resolution. I will use 1 degree of latitude and longitude. (It was originally 2).
res(grid) <- 1

# Make the grid have the same coordinate reference system (CRS) as the shapefile.

# Transform this raster into a polygon to get a grid, but without the shapefile.
gridpolygon <- rasterToPolygons(grid)

# Intersect the grid with the Neotropical shapefile. 
neo.grid <- intersect(neo, gridpolygon)

# Plot the intersected shape to see if everything is fine.

# add the points corresponding to species occurrences
points(alagl$lon, alagl$lat, col='orange', pch=20, cex=0.75)
# plot points again to add a border for better visibility
points(alagl$lon, alagl$lat, col='red', cex=0.75)

Here is the map I got:

enter image description here

And here is an example of map I am looking for (source: DOI: 10.1371/journal.pone.0152468):

enter image description here


1 Answer 1


Here's the approach I would take. The general idea is that we create a grid over your countries, then we do a spatial join to the grid then count the number of observations in each grid.

In this example, I just count the number of points. This does not at all match what you provided as an example image. There may be a field that has a population count that should be used as a weight.


# read in the polygons
polygons <- read_sf("/path/to/shapefile/Lowenberg_Neto_2014.shp")

# get the population data
alag <- dismo::gbif("Allagoptera", "*", geo = TRUE)

# remove missing locations
# cast to sf object
# transform to 4267 (same CRS as polygons)
species <- tibble::as_tibble(alag) |> 
  filter(!is.na(lat), !is.na(lon)) |> 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) |> 
  st_transform(crs = 4267)

# make a fishnet grid over the countries
grd <- st_make_grid(polygons, n = 50)
# visualize the grid

# find which grid points intersect `polygons` (countries) 
# and create an index to subset from
index <- which(lengths(st_intersects(grd, polygons)) > 0)

# subset the grid to make a fishnet
fishnet <- grd[index]

# visualize the fishnet

# cast the fishnet as an sf object, join the species to it
# count the number of observations by grid ID and plot it 
fishnet |>  
  st_as_sf() |> # cast to sf
  mutate(grid_id = row_number()) |> # create unique ID
  st_join(species) |> # join the species dataset
  group_by(grid_id) |> # group by the grid id
  count() |> # count the number of rows
  # fill the plot by the number of points in the grid
  ggplot(aes(fill = n)) + 
  # make kinda pretty 
  geom_sf(lwd = 0.1, color = "white")

enter image description here

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.