1

I am new to spatial data manipulation and now trying to use R to construct a cell-year panel. I wonder if any of you could kindly provide any advice. I have two main datasets, both geo-referenced. One is a conflict dataset from https://www.acleddata.com/download/2909/. The location of conflict is geocoded so I used the coordinates to create a SpatialPointsDataFrame. The other dataset is similar to this one, where I have the coordinates of the aid project location and other information. Therefore now I have two SpatialPointsDataFrames.

What I want to do are:

  1. get a map of Africa and make my datasets SpatialPointsDataFrame

  2. create 100km*100km grid cells covering entire Africa

  3. mask (or any other jargon) the grid and map of Africa to get the cells that are on land, instead of seas

  4. join (or do a sp::over?) the aid projects and conflicts data to every cell, so that I have a dataframe consisting of all cells in Africa and the information of all aid projects and conflicts within each cell. (so some cells may not have any aid project or conflict and some may have many)

  5. the very important thing is, I want to aggregate all data (both aid and conflict) in a cell by year. The final product will be a dataframe that looks like the following: enter image description here

The variables "count aid" and "count conflict" are the numbers of aid projects and conflicts in a given year.

I have tried a lot to construct this cell-year panel. Here are what I have done:

Step 1 Prepare data for manipulation

#load the libraries
library(sf)
library(sp) 
library(raster)
library(dplyr)
library(rworldmap)
library(rmapshaper)
setwd("D:/")

#Step 1 get the map of Africa
isoall <- ccodes()
isoafrica <- filter(isoall, isoall$CONTINENT == "Africa")
isoafrica <- isoafrica[, "ISO3"]

#now loop through African countries and dump each into a list
#I borrowed this code from obrl_soil's answer in https://gis.stackexchange.com/questions/228706/associating-country-to-grid-cell-in-r

GADM0 <- list()
for (c in isoafrica){
  GADM0[[c]] <- getData("GADM", country=c, level=0)
}

africa <- do.call('rbind', GADM0)
object.size(africa)

#simplify the geometry, as obrl_soil suggests 
africa <- ms_simplify(africa, keep = 0.05) 
object.size(afirca)
str(africa, max.level = 2)
class(africa)
africa <- africa[, 1:4]

#project the map of Africa
crs <- "+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
africa <- spTransform(africa, crs)

#I use some sample data to create a SpatialPointsDataFrame
conflict <- data.frame("eventid" = c("1", "2", "7890", "41232", "9803"),
                       "year" = c(2003, 2008, 2013, 2002, 2006),
                       "lat" = c(-14.72278, -16.15639, -12.97395, -17.49194, -25.41444),
                       "long" = c(34.36083, 33.58667, 40.51775, 37.02889, 32.58611),
                       "fatalities" = c(12, 31, 44, 4, 0))

# Add Coordinate Reference System to SPDF objects    
coordinates(conflict) <- c("long", "lat")
wgs84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(conflict) <- wgs84
conflict <- spTransform(conflict, crs)

I have 2 practical questions here:

a. I used this CRS (Africa Albers Equal Area Conic; ESRI:102022) to project the map and other objects ("+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs") is this CRS suitable for my work?

b. When I project 2 SpatialPointsDataFrames, I found I have to first assign an unprojected CRS and then I can spTransform it to a projected CRS. Am I correct?

Step 2 create 100km*100km grid cells covering entire Africa

I tried very hard to find a way to create grid cells. I guess both raster and SpatialGridDataFrame work. But the latter one is more intuitive for me. But later I realized that I may have to convert it to a raster to do a mask. Here is my code.

#create a SpatialGridDataFrame
bb <- bbox(africa)
cs <- c(100000, 100000)  # cell size 
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells in each direction
grid <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
sp_grid <- SpatialGridDataFrame(grid, data=data.frame(id=1:prod(cd)),
                               proj4string=crs)

Step 3: mask the grid and map of Africa to get the cells that are on land

#I failed to do an over of Africa on a SpatialGridDataFrame, so I convert the SPDF to a raster for mask
raster_grid <- raster(extent(sp_grid))
res(raster_grid) <- c(100000, 100000)
proj4string(raster_grid) <- crs
raster_grid[] <- runif(ncell(raster_grid), 1, 10) #assign some random value to the raster
mask_grid <- mask(raster_grid, africa) 
plot(mask_grid)
newgrid <- rasterToPolygons(mask_grid, n = 4, na.rm = TRUE)

From the plot, I think I got it right. Now I have the grid of all cells on the land.

I feel it is especially hard to do step 4 and 5. It is easy to find which aid project/conflict falls in a cell by sp::over, but then it returned a list with a length of the aid/conflict data. what I want is a dataframe with the cell id in the first column, so I tried the following code.

conflictingrid <- sp::over(x = newgrid, y = conflict, returnList = TRUE)

I am not sure if I am right and don't know what to do next.

Can anyone help me in Step 4 and 5?

Your Answer

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

Browse other questions tagged or ask your own question.