barycentric examples

[mdsumner]

These just to keep me on track, very experimental

[mdsumner]

library(raster)
e <-new("Extent", xmin = 481.897772388206, xmax = 934.102227611794, 
         ymin = 47, ymax = 392)
p_srs <- "+proj=stere +lon_0=165 +lat_0=-90 +lat_ts=-71 +datum=WGS84"
polar <- raster(extent(-4e6, 5.4e6, 0, 7e6), crs = p_srs, res = 25000)

library(angstroms)
f <- raadtools::cpolarfiles()$fullname[1]
temp <- crop(romsdata(f, varname = "temp", slice = c(1, 31)), e)
coords <- crop(romscoords(f, transpose = TRUE), e)

index <- bary_index(temp, coords = coords, grid = polar)
ok <- !is.na(index$idx)
r <- setValues(polar, NA_integer_)
r[ok] <- colSums(matrix(values(temp)[index$tri[, index$idx[ok]]], nrow = 3) * t(index$p)[, ok])
plot(r, asp = 1)

[mdsumner]

Note the need to expand out the xt/yt_ocean and flip:

library(raster)
e <-new("Extent", xmin = 0, xmax = 360, 
        ymin = 0, ymax = 150)
p_srs <- "+proj=stere +lon_0=165 +lat_0=-90 +lat_ts=-71 +datum=WGS84"
polar <- raster(extent(-4e6, 5.4e6, 0, 7e6), crs = p_srs, res = 8000)

library(angstroms)
f <- "/perm_storage/data/ocean_bgc.nc-20910630"

x <- crop(romsdata(f, varname = "fe", slice = c(1, 1)), e)
#coords <- crop(romscoords(f,  c("geolon_t", "geolat_t"), transpose = TRUE), e)
xx <- rawdata(f, "xt_ocean")
yy <- rawdata(f, "yt_ocean")[1:150]
coords <- flip(setValues(brick(x, x), as.matrix(expand.grid(xx, yy))), "y")
index <- bary_index(x, coords = coords, grid = polar)
ok <- !is.na(index$idx)
r <- setValues(polar, NA_integer_)
r[ok] <- colSums(matrix(values(x)[index$tri[, index$idx[ok]]], nrow = 3) * t(index$p)[, ok])
plot(r, asp = 1)
plot(spTransform(raadtools::polar_map(), p_srs), add = TRUE)

[mdsumner]

Now a file with full x/y curvilinear (though I think these are just xt/yt as above, need to check.

library(raster)
e <-new("Extent", xmin = 0, xmax = 360, 
        ymin = 0, ymax = 150)
p_srs <- "+proj=stere +lon_0=165 +lat_0=-90 +lat_ts=-71 +datum=WGS84"
polar <- raster(extent(-4e6, 5.4e6, 0, 7e6), crs = p_srs, res = 8000)

library(angstroms)
f <- "/perm_storage/data/ocean_month.nc-20910630"
x <- crop(romsdata(f, varname = "sea_level", slice = c(1, 1)), e)
coords <- crop(romscoords(f,  c("geolon_t", "geolat_t"), transpose = TRUE), e)

index <- bary_index(x, coords = coords, grid = polar)
ok <- !is.na(index$idx)
r <- setValues(polar, NA_integer_)
r[ok] <- colSums(matrix(values(x)[index$tri[, index$idx[ok]]], nrow = 3) * t(index$p)[, ok])
plot(r, asp = 1, col = viridis::viridis(100))
plot(spTransform(raadtools::polar_map(), p_srs), add = TRUE)

ping @SnowpeaSoho FYI this is with remotes::install_github("hypertidy/quadmesh")

[mdsumner]

still WIP but close now

library(raster)
xthing <- readRDS("mike_1.rds")
## coords is longitude,latitude in a 2-layer RasterStack/Brick
coords <- xthing$cc
## x is the data, 'tos' for example
x <- xthing$tos
## two helper functions for raster
index_coords <- function(x) {
cell <- seq_len(raster::ncell(x))
raster::setValues(x, cbind(raster::colFromCell(x, cell),
raster::rowFromCell(x, rev(cell))))
}
index_r <- function(x) {
raster::setExtent(x, raster::extent(0, ncol(x), 0, nrow(x)))
}
## helper fun to wrap longitude into -180,180
make_180 <- function(x) {
a360 <- x[,1] > 180
x[a360,1] <- x[a360,1] - 360
x
}
## helper fun to get longlat from grid
longlatFromGrid <- function(x) {
if (raster::isLonLat(x)) return(raster::xyFromCell(x, seq_len(raster::ncell(x))))
xy <- xyFromCell(x, 1:raster::ncell(x))
as.matrix(reproj::reproj(xy, "OGC:CRS84", source = raster::projection(x))[, 1:2, drop = FALSE])
}
## barycentric lookup in index-space
index <- quadmesh::bary_index(x, coords = index_coords(coords), grid = index_r(x))
ok <- !is.na(index$idx)
grid <- raster::raster(raster::extent(-180, 180, -90, 90), ncols = 360, nrows = 180)
d <- pi * 6378137
grid <- raster(extent(-1, 1, -1/2, 1/2) * d, nrows = 180, ncols = 360, crs = "+proj=moll")
grid <- raster(extent(-1, 1, -1, 1) * 2e7, nrows = 360, ncols = 360, crs = "+proj=omerc +lon_0=147 +gamma=10 +lonc=180")
grid <- raster(extent(-1, 1, -1, 1) * 2e7, nrows = 360, ncols = 360, crs = "+proj=tmerc +lon_0=147")
## map between target grid coordinates and the index grid
XY <- raster::values(coords)
#LL <- xyFromCell(grid, 1:ncell(grid))
LL <- longlatFromGrid(grid)
bad <- !is.finite(LL[,1]) | is.na(LL[,2])
nn <- RANN::nn2(make_180(XY)[ok, ], LL[!bad, ], k = 1)
iidx <- seq_len(ncell(grid))
iidx[bad] <- NA
r <- setValues(grid, NA_integer_)
r[na.omit(iidx)] <- colSums(matrix(values(x)[index$tri[, index$idx[ok]]], nrow = 3) * t(index$p)[, ok])[nn$nn.idx[,1]]
library(raster)
plot(r, col =hcl.colors(25))

https://twitter.com/mdsumner/status/1470285336106323969?s=20

[mdsumner]

and now as functions

library(raster)
xthing <- readRDS("mike_1.rds")

## coords is longitude,latitude in a 2-layer RasterStack/Brick
coords <- xthing$cc
## x is the data, 'tos' for example
x <- xthing$tos

## two helper functions for raster
index_coords <- function(x) {
  cell <- seq_len(raster::ncell(x))
  raster::setValues(x, cbind(raster::colFromCell(x, cell),
                             raster::rowFromCell(x, rev(cell))))
}
index_r <- function(x) {
  raster::setExtent(x, raster::extent(0, ncol(x), 0, nrow(x)))
}
## helper fun to wrap longitude into -180,180
make_180 <- function(x) {
  a360 <- x[,1] > 180
  x[a360,1] <- x[a360,1] - 360
  x
}
## helper fun to get longlat from grid
longlatFromGrid <- function(x) {
  if (raster::isLonLat(x)) return(raster::xyFromCell(x, seq_len(raster::ncell(x))))
  xy <- xyFromCell(x, 1:raster::ncell(x))
  suppressWarnings(rgdal::project(xy, raster::projection(x), inv = TRUE))
}

## this is our target grid, pick regular longlat or whatever you want
grid <- raster::raster(raster::extent(-180, 180, -90, 90), ncols = 360, nrows = 180)
d <- pi * 6378137
grid <- raster(extent(-1, 1, -1/2, 1/2) * d, nrows = 180, ncols = 360, crs = "+proj=moll")
#grid <- raster(extent(-1, 1, -1, 1) * 2e7, nrows = 360, ncols = 360, crs = "+proj=omerc +lon_0=147 +gamma=10 +lonc=180")
#grid <- raster(extent(-1, 1, -1, 1) * 2e7, nrows = 360, ncols = 360, crs = "+proj=tmerc +lon_0=147")

raster0 <- function(x, ...) {
  x <- raster::raster(x, ..., stopIfNotEqualSpaced = FALSE)
  raster::setExtent(x, raster::extent(0, ncol(x), 0, nrow(x)))
}
regrid <- function(x, coords = NULL, grid = NULL,  ...) {
  UseMethod("regrid")
}

regrid.character <- function(x, coords = NULL, grid = NULL,  ...) {
  if (!is.null(coords)) {
    if (inherits(coords, "BasicRaster")) {
      #coords <- raster::brick(raster0(x, varname = coords[1]), raster0(x, coords[2]))
    } else {
      coords <- raster::brick(raster0(x, varname = coords[1]), raster0(x, varname = coords[2]))
    }
  }
  ## here try different names
  if (is.null(coords)) {
    coords <- raster::brick(raster0(x, varname = "lon"), raster0(x, varname = "lat"))
  }
  regrid(raster0(x,...), coords = coords, grid = grid)
}
regrid.BasicRaster <- function(x, coords = NULL, grid = NULL, ...) {
  ## barycentric lookup in index-space
  if (is.null(coords)) stop("coords must be set for raster input")
  index <- quadmesh::bary_index(x, coords = index_coords(coords), grid = index_r(x))
  ok <- !is.na(index$idx)
  ## map between target grid coordinates and the index grid
  XY <- raster::values(coords)

  if (is.null(grid)) grid <- raster::raster(ncols = 360, nrows = 180)
  LL <- longlatFromGrid(grid)
  bad <- !is.finite(LL[,1]) | is.na(LL[,2])
  nn <- RANN::nn2(make_180(XY)[ok, ], LL[!bad, ], k = 1)

  iidx <- seq_len(ncell(grid))
  iidx[bad] <- NA
  r <- setValues(grid, NA_integer_)
  r[iidx[!is.na(iidx)]] <- colSums(matrix(values(x)[index$tri[, index$idx[ok]]], nrow = 3) * t(index$p)[, ok])[nn$nn.idx[,1]]
 r
}