gdalio
gdalio copied to clipboard
do this helper stuff for SDS for mesh plots
rinfo <- function(x, sds = NULL) vapour::vapour_raster_info(x, sds = sds)
size <- c(512, 512)
#' Read quickly via gdal
#'
#' @param x
#' @param size
#' @param sds
#'
#' @return
#' @export
#'
#' @examples
gdal_approx <- function(x, size = c(128, 128), sds = NULL) {
vi <- rinfo(x, sds = sds)
if (vi$extent[3] > vi$extent[4]) vi$extent[3:4] <- vi$extent[4:3]
if (is.null(size)) size <- vi$dimXY
val <- vapour::vapour_read_raster_dbl(x, sds = sds, window = c(0,0, vi$dimXY, size))
raster::raster(t(matrix(val, size[1])[,size[2]:1]), xmn = vi$extent[1], xmx = vi$extent[2],
ymn = vi$extent[3], ymx = vi$extent[4])
}
fp <- function(x) {
file.path("S3A_OL_2_WFR/S3A_OL_2_WFR____20220401T141026_20220401T141326_20220403T010132_0179_083_338_3600_MAR_O_NT_003.SEN3",
x)
}
file <- fp("Oa01_reflectance.nc")
xyfile <- fp("geo_coordinates.nc")
.data <- gdal_approx(file, sds = 1)
plot(.data)
.xydata <- brick(gdal_approx(xyfile, sds = 3),
gdal_approx(xyfile, sds = 2)) * 1.e-06
library(anglr)
mesh_plot(.data, coords = .xydata)
maps::map(add = TRUE)
axis(1); axis(2)
also see https://github.com/hypertidy/vapour/issues/118