'union' between spatial polygons and spatial points in R

Question

I have a gridded map shapefile with grid identifier cells (gid).

dput(head(grid))
structure(list(gid = c("183253", "183254", "183255", "183256", 
"183258", "183259"), xcoord = c(6.25, 6.75, 7.25, 7.75, 8.75, 
9.25), ycoord = c(37.25, 37.25, 37.25, 37.25, 37.25, 37.25), 
    col = c("373", "374", "375", "376", "378", "379"), row = c("255", 
    "255", "255", "255", "255", "255"), gwno = c(615L, 615L, 
    615L, 615L, 616L, 616L), country = c("Algeria", "Algeria", 
    "Algeria", "Algeria", "Tunisia", "Tunisia"), km2 = c(3875.43674549, 
    3875.43674549, 3875.43674549, 3875.43674549, 3875.43674549, 
    3875.43674549), geometry = structure(list(structure(list(
        structure(c(6, 6, 6.5, 6.5, 6, 37, 37.5, 37.5, 37, 37
        ), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg"
    )), structure(list(structure(c(6.5, 6.5, 7, 7, 6.5, 37, 37.5, 
    37.5, 37, 37), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", 
    "sfg")), structure(list(structure(c(7, 7, 7.5, 7.5, 7, 37, 
    37.5, 37.5, 37, 37), .Dim = c(5L, 2L))), class = c("XY", 
    "POLYGON", "sfg")), structure(list(structure(c(7.5, 7.5, 
    8, 8, 7.5, 37, 37.5, 37.5, 37, 37), .Dim = c(5L, 2L))), class = c("XY", 
    "POLYGON", "sfg")), structure(list(structure(c(8.5, 8.5, 
    9, 9, 8.5, 37, 37.5, 37.5, 37, 37), .Dim = c(5L, 2L))), class = c("XY", 
    "POLYGON", "sfg")), structure(list(structure(c(9, 9, 9.5, 
    9.5, 9, 37, 37.5, 37.5, 37, 37), .Dim = c(5L, 2L))), class = c("XY", 
    "POLYGON", "sfg"))), class = c("sfc_POLYGON", "sfc"), precision = 0, bbox = structure(c(xmin = 6, 
    ymin = 37, xmax = 9.5, ymax = 37.5), class = "bbox"), crs = structure(list(
        input = "WGS 84", wkt = "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"latitude\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"longitude\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(gid = NA_integer_, 
xcoord = NA_integer_, ycoord = NA_integer_, col = NA_integer_, 
row = NA_integer_, gwno = NA_integer_, country = NA_integer_, 
km2 = NA_integer_), .Label = c("constant", "aggregate", "identity"
), class = "factor"), row.names = c(NA, 6L), class = c("sf", 
"data.frame"))

I have spatial points that represent the sum of fatalities during a period of time.

 dput(head(acled))
structure(list(OBJECTID = c(1, 2, 3, 4, 5, 6), data_id = c(8037001, 
6225870, 8037590, 8037835, 6270986, 8037919), year = c(2010, 
2010, 2009, 2007, 2007, 2007), event_type = c("Battles", "Battles", 
"Battles", "Battles", "Battles", "Battles"), country = c("Angola", 
"Angola", "Angola", "Angola", "Angola", "Angola"), latitude = c(-8.8383, 
-5.5758, -4.7666, -4.7666, -5.5758, -4.7666), longitude = c(13.2344, 
12.1871, 12.55, 12.55, 12.1871, 12.55), fatalities = c(2, 0, 
0, 1, 2, 1), geometry = structure(list(structure(c(13.2344, -8.8383
), class = c("XY", "POINT", "sfg")), structure(c(12.1871, -5.5758
), class = c("XY", "POINT", "sfg")), structure(c(12.55, -4.7666
), class = c("XY", "POINT", "sfg")), structure(c(12.55, -4.7666
), class = c("XY", "POINT", "sfg")), structure(c(12.1871, -5.5758
), class = c("XY", "POINT", "sfg")), structure(c(12.55, -4.7666
), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = 12.1871, ymin = -8.8383, 
xmax = 13.2344, ymax = -4.7666), class = "bbox"), crs = structure(list(
    input = "WGS 84", wkt = "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"latitude\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"longitude\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(OBJECTID = NA_integer_, 
data_id = NA_integer_, year = NA_integer_, event_type = NA_integer_, 
country = NA_integer_, latitude = NA_integer_, longitude = NA_integer_, 
fatalities = NA_integer_), .Label = c("constant", "aggregate", 
"identity"), class = "factor"), row.names = c(NA, 6L), class = c("sf", 
"data.frame"))

Ideally, I would like a final spatialpolygonsdataframe that represents polygons summing the fatalities inside them for each year + the polygon that don't have any fatalities. It would be like the 'union' function between two similar shapefiles.

However, when I try do so with the function point.in.poly(acled, grid) from the spatialEco package I get a spatialpointsdataframe with only the points falling into the corresponding grid polygons. Same with the 'over' function.

Is there a way to get the opposite ? I tried over(grid, ACLED) but it does not give me a polygon per year as there is a total number of 10667 cells in the grid shapefile, in other words, it does not duplicate the polygons to let me distinguish the fatalities for the different years.

Answer 1

I found out how to do it.

Two ways:

1. applying point.in.poly(grid, acled). So x= polygons, y=points.

2. applying point.in.poly(ACLED, grid) in the right way then merging with the initial grid.

grid.acled <- sp::merge(grid, pts.poly[c("fatalities","event_type", "gid", "year")], by="gid", all.x= TRUE).