1

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
library(rgdal)
library(raster)
library(rgeos)
library(dismo)

#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:
#https://rfunctions.blogspot.com/2014/12/how-to-create-grid-and-intersect-it.html
#I just replaced the shapefile with the one for the Neotropics and modified the comments accordingly.


Load the Neotropical shapefile.


#Data come from here:
#https://figshare.com/articles/dataset/Neotropical_region_a_shapefile_of_Morrone_s_2014_biogeographical_regionalisation/3569361


#Function 'download.shapefile' comes from here: 
#https://www.r-bloggers.com/2013/09/an-r-function-to-download-shapefiles/


download.shapefile("https://figshare.com/articles/dataset/Neotropical_region_a_shapefile_of_Morrone_s_2014_biogeographical_regionalisation/3569361?file=5646714",
                          "Lowenberg_Neto_2014")


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


Plot the shape to see if everything is fine.


plot(neo)


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.


proj4string(grid)<-proj4string(neo)


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.


plot(neo.grid)


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

Jhonny
  • 11
  • 1
  • 3

1 Answers1

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.

library(sf)
library(dplyr)
library(ggplot2)

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


plot(grd)


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


plot(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

thus__
  • 179
  • 6