reprojection scheme for vector field data

Question

reprojection scheme for vector field data

Closed this issue 5 years ago · 5 comments

This seems to work for converting surface currents to a polar projection.

The key is to create "end points" from the du/dv values in Plate Carree (equirectangular), then transform the start and end points into polar, then recreate the du/dv values, and use nearest neighbour to rebuild a grid. Each step can be indexed for pretty fast implementation.

library(raadtools)
library(dplyr)
## U/V components in longlat
u <- readcurr(xylim = extent(-180, 180, -75, -20), uonly = TRUE)
v <- readcurr(xylim = extent(-180, 180, -75, -20), vonly = TRUE)

## target grid, 24km polar stereographi
target <- sospatial:::stere_grid(24)
library(tabularaster)
## get the values and cell index on a data frame
d <- as_tibble(u) %>% rename(u = cellvalue)
d$v <- values(v)

## project to equrectangular, a projection orthogonal to the U/V vector components
d[c("x", "y")] <- rgdal::project(xyFromCell(u, d$cellindex), "+proj=eqc")
## update x and y by u and v (could be scaled if needed)
d$x1 <- d$x + d$u
d$y1 <- d$y +  d$v

## project the starting points from the original longlat to the target
d[c("sx", "sy")] <- rgdal::project(xyFromCell(u, d$cellindex), projection(target))
## unproject the end points from eqc to longlat and then forward to the target
d[c("ex", "ey")] <- rgdal::project(rgdal::project(as.matrix(d[c("x1", "y1")]), "+proj=eqc", inv = TRUE), 
                                   projection(target))

## update to projected u and projected v
d$pu <- d$sx - d$ex
d$pv <- d$sy - d$ey

## regrid by nn
library(nabor)
xy <- coordinates(target)
knn <- WKNNF(as.matrix(d[c("sx", "sy")]))
idx <- knn$query(xy, k = 1, eps = 0)
uu <- setValues(target, d$pu[idx$nn.idx])
vv <- setValues(target, d$pv[idx$nn.idx])
plot(sqrt(uu * uu + vv * vv))

target can be made arbitrarily

TODO: need to re-mask land and areas north of the maximum (or trim the query coordinates)

Answer 1 · 2018-05-08T07:31:18.000Z

This is better, no need for masking after the fact

library(raadtools)
library(dplyr)
u <- readcurr(xylim = extent(-180, 180, -75, -20), uonly = TRUE)
v <- readcurr(xylim = extent(-180, 180, -75, -20), vonly = TRUE)

target <- sospatial:::stere_grid(24)
library(tabularaster)
d <- as_tibble(u) %>% rename(u = cellvalue)
d$v <- values(v)
d[c("x", "y")] <- rgdal::project(xyFromCell(u, d$cellindex), "+proj=eqc")
d$x1 <- d$x + d$u
d$y1 <- d$y +  d$v

d[c("sx", "sy")] <- rgdal::project(xyFromCell(u, d$cellindex), projection(target))
d[c("ex", "ey")] <- rgdal::project(rgdal::project(as.matrix(d[c("x1", "y1")]), "+proj=eqc", inv = TRUE), 
                                   projection(target))


d$pu <- d$sx - d$ex
d$pv <- d$sy - d$ey
dim(d)



library(nabor)
xy <- coordinates(target)
ee <- extract(u, rgdal::project(xy, projection(target), inv = TRUE))

xy <- xy[!is.na(c(ee)), ]
knn <- WKNNF(as.matrix(d[c("sx", "sy")]))
idx <- knn$query(xy, k = 1, eps = 0)
uu <- vv <- target
uu[!is.na(ee)] <- d$pu[idx$nn.idx]
vv[!is.na(ee)] <- d$pv[idx$nn.idx]
plot(sqrt(uu * uu + vv * vv), col = palr::sstPal(64)[1:40])

Answer 2 · 2018-05-08T15:11:23.000Z

Here's the efficient-ish job

default_grid <- function() {
  prjj <-         "+proj=laea +lat_0=-90 +datum=WGS84"
  raster(spex::buffer_extent(projectExtent(raster(extent(-180, 180, -90, -30), 
                                                  crs = "+init=epsg:4326"), 
                                           prjj), 25000), 
         res = 25000, crs = prjj)
  
  
}


library(raadtools)
library(dplyr)
files <- currentsfiles()

#u <- readcurr(xylim = extent(-180, 180, -75, -20), uonly = TRUE)
#v <- readcurr(xylim = extent(-180, 180, -75, -20), vonly = TRUE)

ex <- extent(-180, 180, -75, -30)
geti_u <- function(i = 1) {
  readcurr(files$date[i], xylim = ex, uonly = TRUE, inputfiles = files)
}
geti_v <- function(i = 1) {
  readcurr(files$date[i], xylim = ex, vonly = TRUE, inputfiles = files)
}
u <- geti_u()
target <- default_grid()
library(nabor)
xy <- coordinates(target)
xyi <- rgdal::project(xy, projection(target), inv = TRUE)
library(tabularaster)
index <- tibble(cellindex = seq_len(ncell(u)))
index[c("qx", "qy")] <- rgdal::project(xyFromCell(u, index$cellindex), "+proj=eqc")
index[c("sx", "sy")] <- rgdal::project(xyFromCell(u, index$cellindex), projection(target))
knn <- WKNNF(as.matrix(index[c("sx", "sy")]))

i <- 1

for (i in seq_along(lv)) {
U <- geti_u(i)
V <- geti_v(i)
index$x1 <- index$qx + values(U)
index$y1 <- index$qy + values(V)

index[c("ex", "ey")] <- rgdal::project(rgdal::project(as.matrix(index[c("x1", "y1")]), "+proj=eqc", inv = TRUE), 
                                   projection(target))

index$pu <- index$sx - index$ex
index$pv <- index$sy - index$ey

ee <- extract(U, xyi)

xyq <- xy[!is.na(c(ee)), ]
idx <- knn$query(xyq, k = 1, eps = 0)
uu <- vv <- target
uu[!is.na(ee)] <- index$pu[idx$nn.idx]
vv[!is.na(ee)] <- index$pv[idx$nn.idx]
writeRaster(uu, filename = sprintf("data-local/aad.gov.au/currents/fileu_%05i.grd", i), overwrite = TRUE)
writeRaster(vv, filename = sprintf("data-local/aad.gov.au/currents/filev_%05i.grd", i), overwrite = TRUE)

}

Bundle them into a VRT (and so forth) with

gdalbuildvrt u.vrt -separate data-local/currents/fileu*.grd

Answer 3 · 2018-09-12T01:42:02.000Z

Another update, better naming and faster to update

  prjj <-         "+proj=laea +lat_0=-90 +datum=WGS84"
  raster(spex::buffer_extent(projectExtent(raster(extent(-180, 180, -90, -30), 
                                                  crs = "+init=epsg:4326"), 
                                           prjj), 25000), 
         res = 25000, crs = prjj)
  
  
}


library(raadtools)
library(dplyr)
files <- currentsfiles()

#u <- readcurr(xylim = extent(-180, 180, -75, -20), uonly = TRUE)
#v <- readcurr(xylim = extent(-180, 180, -75, -20), vonly = TRUE)

ex <- extent(-180, 180, -75, -30)
geti_u <- function(i = 1) {
  readcurr(files$date[i], xylim = ex, uonly = TRUE, inputfiles = files)
}
geti_v <- function(i = 1) {
  readcurr(files$date[i], xylim = ex, vonly = TRUE, inputfiles = files)
}
u <- geti_u()
target <- default_grid()
library(nabor)
xy <- coordinates(target)
xyi <- rgdal::project(xy, projection(target), inv = TRUE)
library(tabularaster)
index <- tibble(cellindex = seq_len(ncell(u)))
index[c("qx", "qy")] <- rgdal::project(xyFromCell(u, index$cellindex), "+proj=eqc")
index[c("sx", "sy")] <- rgdal::project(xyFromCell(u, index$cellindex), projection(target))
knn <- WKNNF(as.matrix(index[c("sx", "sy")]))

i <- 1

#gsub( "nc$", "grd", gsub("dt_global_allsat_phy_l4_", "dt_south_polar_u_", basename(files$fullname[i])))
dirbase <- "/rdsi/PRIVATE/raad/data_local/aad.gov.au/currents/polar"
for (i in seq_len(nrow(files))) {
  ufile <- gsub( "nc$", "grd", gsub("dt_global_allsat_phy_l4_", "dt_south_polar_u_", basename(files$fullname[i])))
  vfile <- gsub( "nc$", "grd", gsub("dt_global_allsat_phy_l4_", "dt_south_polar_v_", basename(files$fullname[i])))
 if (file.exists(ufile) && file.exists(vfile)) next; 
  U <- geti_u(i)
  V <- geti_v(i)
  index$x1 <- index$qx + values(U)
  index$y1 <- index$qy + values(V)
  
  index[c("ex", "ey")] <- rgdal::project(rgdal::project(as.matrix(index[c("x1", "y1")]), "+proj=eqc", inv = TRUE), 
                                         projection(target))
  
  index$pu <- index$sx - index$ex
  index$pv <- index$sy - index$ey
  
  ee <- extract(U, xyi)
  
  xyq <- xy[!is.na(c(ee)), ]
  idx <- knn$query(xyq, k = 1, eps = 0, radius = 0)
  uu <- vv <- target
  uu[!is.na(ee)] <- index$pu[idx$nn.idx]
  vv[!is.na(ee)] <- index$pv[idx$nn.idx]
  ufile <- gsub( "nc$", "grd", gsub("dt_global_allsat_phy_l4_", "dt_south_polar_u_", basename(files$fullname[i])))
  vfile <- gsub( "nc$", "grd", gsub("dt_global_allsat_phy_l4_", "dt_south_polar_v_", basename(files$fullname[i])))
  uu <- setZ(uu, files$date[i])
  vv <- setZ(vv, files$date[i])
  tmp <- writeRaster(uu, filename = file.path( dirbase, ufile), overwrite = TRUE)
  tmp <- writeRaster(vv, filename = file.path( dirbase, vfile), overwrite = TRUE)
  print(i)
}

Answer 4 · 2018-11-05T22:14:45.000Z

final_final_final_really_meanit.R

default_grid <- function() {
  
  prjj <-         "+proj=laea +lat_0=-90 +datum=WGS84"
  raster(spex::buffer_extent(projectExtent(raster(extent(-180, 180, -90, -30), 
                                                  crs = "+init=epsg:4326"), 
                                           prjj), 25000), 
         res = 25000, crs = prjj)
  
  
}


library(raadtools)
library(dplyr)
files <- currentsfiles()

#u <- readcurr(xylim = extent(-180, 180, -75, -20), uonly = TRUE)
#v <- readcurr(xylim = extent(-180, 180, -75, -20), vonly = TRUE)

ex <- extent(-180, 180, -75, -30)
geti_u <- function(i = 1) {
  readcurr(files$date[i], xylim = ex, uonly = TRUE, inputfiles = files)
}
geti_v <- function(i = 1) {
  readcurr(files$date[i], xylim = ex, vonly = TRUE, inputfiles = files)
}
u <- geti_u()
target <- default_grid()
library(nabor)
xy <- coordinates(target)
xyi <- rgdal::project(xy, projection(target), inv = TRUE)
library(tabularaster)
index <- tibble(cellindex = seq_len(ncell(u)))
index[c("qx", "qy")] <- rgdal::project(xyFromCell(u, index$cellindex), "+proj=eqc")
index[c("sx", "sy")] <- rgdal::project(xyFromCell(u, index$cellindex), projection(target))
knn <- WKNNF(as.matrix(index[c("sx", "sy")]))

dirbase <- "/rdsi/PRIVATE/raad/data_local/aad.gov.au/currents/polar"
#for (i in seq_len(nrow(files))) {
for (i in seq_len(nrow(files))) {
  ## split by year
  diryear <- file.path(dirbase, format(files$date[i], "%Y"))
  if (!file.exists(diryear)) dir.create(diryear)
  ufile <- file.path(diryear, gsub( "nc$", "grd", gsub("dt_global_allsat_phy_l4_", "dt_south_polar_u_", basename(files$fullname[i]))))
  vfile <- file.path(diryear, gsub( "nc$", "grd", gsub("dt_global_allsat_phy_l4_", "dt_south_polar_v_", basename(files$fullname[i]))))
  if (file.exists(ufile) && file.exists(vfile)) next; 
  U <- geti_u(i)
  V <- geti_v(i)
  index$x1 <- index$qx + values(U)
  index$y1 <- index$qy + values(V)
  
  index[c("ex", "ey")] <- rgdal::project(rgdal::project(as.matrix(index[c("x1", "y1")]), "+proj=eqc", inv = TRUE), 
                                         projection(target))
  
  index$pu <- index$sx - index$ex
  index$pv <- index$sy - index$ey
  
  ee <- extract(U, xyi)
  
  xyq <- xy[!is.na(c(ee)), ]
  idx <- knn$query(xyq, k = 1, eps = 0, radius = 0)
  uu <- vv <- target
  uu[!is.na(ee)] <- index$pu[idx$nn.idx]
  vv[!is.na(ee)] <- index$pv[idx$nn.idx]
  uu <- setZ(uu, files$date[i])
  vv <- setZ(vv, files$date[i])
  tmp <- writeRaster(uu, filename = ufile, overwrite = TRUE)
  tmp <- writeRaster(vv, filename = vfile, overwrite = TRUE)
  print(i)
}

Answer 5 · 2019-11-25T01:14:02.000Z

This is now "polar_currents/" in https://github.com/AustralianAntarcticDivision/raad-deriv