Aggregating point data regarding lat/lon and year

Question

I have two datasets:

1. The first one contains georeferenced point data. The structure looks like that:

2. The second one contains georeferenced point data + radius. The structure looks like that:

My goal is to get the sum of the "newusdco" for the points in the corresponding area of the second dataset.

My first problem is, that I don't know how to aggregate the data to the lat/long+radius data.

My second problem is, however, that the data varies over time too. Meaning that I only want to get the sum of "newusdco" in the conflict site in the same years.

I have been reading up in the forum a lot but couldn't find a solution, since I am not very familiar with R.

Answer 1

You can follow the next steps to achieve your goal:

1. Transform your data as a data.frame object
2. Convert the data.frame to SpatialPointsDataFrame (SPDF) object
3. Add a Coordinate Reference System (CRS) to SPDF objects
4. Build a circular buffer around the points of the second dataset with the given radius
5. Query which points of the first dataset are inside the buffer polygon
6. Subset and sum of data by year

Try the commented R code below.

• Goal: get the sum of the "newusdco" for the points in the corresponding area of the second dataset in the same years.

Step 1:

### Load libraries

library('sp') # for handling spatial objects
library('rgeos') # for buffers geometry
library('mapview') # for interactive map viewing in R (so nice!)

### Example data

# Data from question
dataset1 <- data.frame("eventid" = c("112130.01", "112087.01", "115838.01", "115839.01", "115839.02"),
                       "year" = c(1989, 1989, 1989, 1989, 1989),
                       "lat" = c(-8.838333, -15.333333, -8.833333, -15.332072, -16.510087),
                       "long" = c(13.234444, 15, 14.5, 12.660025, 15.435258),
                       "newusdco" = c(201289, 1143688, 1276358, 230453.667, 230453.667))

# Data from question
dataset2 <- data.frame("id" = c("70", "70", "70"),
                       "year" = c(1989, 1990, 1991),
                       "latitude" = c(12.19, 12.10, 11.93),
                       "longitude" = c(38.78, 39.03, 38.04),
                       "radius" = c(400, 300, 400),
                       "conflictarea" = c(455918, 280999, 444130),
                       "conflictsize" = c(530, 530, 530),
                       "conflictterritory" = c("Ethiopia", "Ethiopia", "Ethiopia"))

# Data not from question: created data with the aim of give an example to match points inside buffer polygons of dataset2
dataset1b <- data.frame("eventid" = c("112130.01", "112087.01", "115838.01", "115839.01", "115839.02", "112130.01", "115839.01"),
                       "year" = c(1989, 1989, 1989, 1989, 1989, 1990, 1990),
                       "lat" = c(12.191, 12.101, 11.931, 12.189, 12.990, 12.1911, 12.1891),
                       "long" = c(38.781, 39.031, 38.041, 38.779, 39.029, 38.781, 38.7791),
                       "newusdco" = c(201289, 1143688, 1276358, 230453.667, 230453.667, 201289, 230453.667))

Step 2:

### Convert to SpatialPointsDataFrame (SPDF) object

# Establish coordinates
coordinates(dataset1) <- c("long", "lat")
coordinates(dataset2) <- c("longitude", "latitude")
coordinates(dataset1b) <- c("long", "lat")

Step 3:

# Add Coordinate Reference System to SPDF objects
wgs84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

proj4string(dataset1) <- CRS(wgs84)
proj4string(dataset2) <- CRS(wgs84)
proj4string(dataset1b) <- CRS(wgs84)

### Plot data

spplot(dataset1, zcol = "newusdco", main = "dataset1 plot")
spplot(dataset2, zcol = "radius", main = "dataset2 plot")

# Interactive map
mapview(dataset1, color = "#CC0000") + mapview(dataset2, color = "#73D216") + mapview(dataset1b, color = "#75507B")

- Answer to the first problem "aggregate the data to the lat/long+radius data"

Step 4:

# Find UTM projection for SpatialObject
# Reference: https://stackoverflow.com/questions/9186496/determining-utm-zone-to-convert-from-longitude-latitude
# Also check: https://epsg.io/

long2UTM <- function(long) {
  (floor((long + 180)/6) %% 60) + 1
}

long2UTM(coordinates(dataset2)[,"longitude"]) # UTM zone 37

# UTM
projUTM2 <- "+proj=utm +zone=37 +north +datum=WGS84 +units=m +no_defs"
dataset2UTM <- spTransform(dataset2, CRS(projUTM2))

# Make polygons from dataset2 circle buffer of radius = radius (in meters)
dataset2Buffer <- gBuffer(spgeom = dataset2UTM, byid = TRUE, width = dataset2UTM$radius)

# plot
mapview(dataset1, color = "#CC0000") + mapview(dataset2, color = "#73D216") + mapview(dataset1b, color = "#75507B") + mapview(dataset2Buffer, color = "#8AE234")

Step 5:

# Spatial Query: points in polygon
# Check which points of dataset1 are over dataset2

dataset2Buffer <- spTransform(dataset2Buffer, CRS(wgs84)) # Transform to WGS84

pointsInPolygons1 <- sp::over(x = dataset1, y = dataset2Buffer, returnList = TRUE) # Empty: no points over buffer
pointsInPolygons1b <- sp::over(x = dataset2Buffer, y = dataset1b, returnList = TRUE) # 6 points over buffer

$`1`
eventid year newusdco
1 112130.01 1989 201289.0
4 115839.01 1989 230453.7
6 112130.01 1990 201289.0
7 115839.01 1990 230453.7

$`2`
eventid year newusdco
2 112087.01 1989  1143688

$`3`
eventid year newusdco
3 115838.01 1989  1276358

- Answer to the second problem: "get the sum of 'newusdco' in the conflict site in the same years"

Step 6:

results <- lapply(X = pointsInPolygons1b, FUN = function(x) aggregate(newusdco ~ year, FUN = sum, data = x)
)

    $`1`
  year newusdco
1 1989 431742.7
2 1990 431742.7

$`2`
  year newusdco
1 1989  1143688

$`3`
  year newusdco
1 1989  1276358