distance to ice with terra

[mdsumner]

this is about 10000x faster

disti <- function(date, lon, lat) {
  files <- raadfiles::amsr2_3k_daily_files()
  files <- raadtools:::.processFiles(timedateFrom(date), files, "daily")
  
  ice <- rast(files$fullname[1])
  cl <- disagg(as.contour(ice, levels = 15))
  
 ## get the longer lines, it's not as clean as with nsidc
  gi <- perim(cl) > 5e6
  cl <- cl[gi, ]
  dists <- terra::distance(union(cl), project(vect(cbind(lon, lat), crs = "OGC:CRS84"), crs(cl)))
  as.vector(apply(dists, 2, min))
}

dist_ice <- function(xyt) {
  xyt$g <- cut(xyt[[3]], "1 day")
  l <- split(xyt, xyt$g)[unique(xyt$g)]
  for (i in seq_along(l)) {
   
   l[[i]]$dist_ice <- disti(l[[i]][[3]][1], l[[i]][[1]], l[[i]][[2]])
  }
  do.call(rbind, l)$dist_ice
}

xyt <-data.frame(147, -38:-72, as.Date("2015-01-02") + rep(1:5, each = 7))
xyt$dist_ice <- dist_ice(xyt)

[mdsumner]

compare tempo, 3k (x) vs 25k (y), differences in the range 29-38km

[mdsumner]

new draft, worked up with Laura

dist_ice1 <- function(xyt, icereader = read_amsr2_3k_ice, levels = 15) {
  #data(lon = , lat = , dat)
  xyt$g <- cut(xyt[[3]], "1 day")
  xyt[["original_order124566"]]  <- 1:nrow(xyt)
  l <- split(xyt, xyt$g)[unique(xyt$g)]
  out <- vector("list", length(l))
  for (i in seq_along(l)) {
    ice <- terra::rast(icereader(l[[i]][[3]][1]))
    cl <- disagg(as.contour(ice, levels = levels))
    
    clp <- as.points(cl)
    v <- vect(cbind(l[[i]][[1]], l[[i]][[2]]), crs = "OGC:CRS84")
    out[[i]] <- as.data.frame(nearest(v, project(clp, "OGC:CRS84")))
    out[[i]]$date <- l[[i]][[3]][1]
    out[[i]]$iceconc <- extract(ice, project(v, crs(ice)))[,2L]
    
  }
  out <- do.call(rbind, out)
  out <- out[order(out[["original_order124566"]]), ]
  out[["original_order124566"]] <- NULL
  out[["from_id"]] <- NULL
  out
}

[mdsumner]

new version with progress and ordering

dist_ice1 <- function(xyt, icereader = read_amsr2_3k_ice, levels = 15) {
  #data(lon = , lat = , dat)
  xyt$g <- cut(xyt[[3]], "1 day")
  xyt[["original_order"]]  <- 1:nrow(xyt)
  l <- split(xyt, xyt$g)[unique(xyt$g)]
  out <- vector("list", length(l))
  pb <- progress::progress_bar$new(total =length(out))
  for (i in seq_along(l)) {
    ice <- terra::rast(icereader(l[[i]][[3]][1]))
    cl <- disagg(as.contour(ice, levels = levels))
    
    clp <- as.points(cl)
    v <- vect(cbind(l[[i]][[1]], l[[i]][[2]]), crs = "OGC:CRS84")
    out[[i]] <- as.data.frame(nearest(v, project(clp, "OGC:CRS84")))
    out[[i]]$date <- l[[i]][[3]][1]
    out[[i]]$iceconc <- extract(ice, project(v, crs(ice)))[,2L]
    out[[i]]$order <- l[[i]][["original_order"]]
    pb$tick()
  }
  out <- do.call(rbind, out)
  out <- out[order(out[["order"]]), ]
  out[["order"]] <- NULL
  out[["from_id"]] <- NULL
  out[["to_id"]] <- NULL
  
  out
}