CODEWITHSUNDEEP
rrasterize
akamini
akamini

Reputation: 13

Rasterize overlapping polygons summarising time intervals

I would like to rasterize overlapping polygons, with the function checking the start/end times for each overlapping polygon and combine them to calculate the total duration of all time periods for those polygons. For example consider 3 overlapping polygons:

Polygon 1: 0800 - 1000 (2 hours) Polygon 2: 0800 - 1200 (4 hours) Polygon 3: 0700 - 0900 (2 hours)

On rasterize of these polygons the output would be 5 hours (0700 - 1200), so the pixels at this location would have a z value of 5.

library(sf)
library(tidyverse)

df <- structure(list(lat = c(53.218461547511, 53.2354543814792, 53.2219368310548,
                             53.2142084420887, 53.2051969857393), lon = c(6.57022923569802,
                                                                          6.55005158685787, 6.57015231400986, 6.5633788874466, 6.5692156138519
                             ), startUTC=c("2023-02-21 07:14:38", "2023-02-21 07:05:37", "2023-02-21 07:54:27", "2023-02-21 17:22:04", "2023-02-21 07:46:53"), endUTC=c("2023-02-21 17:32:16", "2023-02-21 18:07:44", "2023-02-21 08:28:18", "2023-02-21 17:52:57", "2023-02-21 08:17:22")), row.names = c(NA, -5L), class = c("tbl_df", "tbl", "data.frame"
                             ))

circle_df <- df %>%
    mutate(startUTC = ymd_hms(startUTC), endUTC = ymd_hms(endUTC)) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  st_buffer(dist = units::set_units(1, "kilometers"))

Here is some code to simulate 5 overlapping polygons and date/time values from my dataset as startUTC/endUTC values. I have investigated terra::rasterize and raster::rasterize, but I cannot see a way to reconcile the intervals within the allowable function.

I also tried to intersect the polygons using st_intersect(), thinking I could work with the intersected parts and their origins attribute, then rasterize the intersected polygons. However the st_intersect fails with a GEOS error which seems to be known and described here:

text

Upvotes: 0

Views: 46

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47546

With these data

df <- data.frame(
    lat = c(53.218, 53.235, 53.222, 53.214, 53.205), 
    lon = c(6.57, 6.55, 6.57, 6.563, 6.569), 
    startUTC = c("2023-02-21 07:14:38", "2023-02-21 07:05:37", "2023-02-21 07:54:27", "2023-02-21 17:22:04", "2023-02-21 07:46:53"), 
    endUTC = c("2023-02-21 17:32:16", "2023-02-21 18:07:44", "2023-02-21 08:28:18", "2023-02-21 17:52:57", "2023-02-21 08:17:22"
))

You can do

df$hours <- as.numeric(as.POSIXct(df$endUTC) - as.POSIXct(df$startUTC)) / 60

library(terra)
v <- vect(df, crs="+proj=longlat")
b <- buffer(v, 1000)

x <- rast(b, res=.001)
x <- rasterize(b, x, "hours", fun="sum")

Upvotes: 0

Related Questions