gdalio icon indicating copy to clipboard operation
gdalio copied to clipboard

do this helper stuff for SDS for mesh plots

Open mdsumner opened this issue 3 years ago • 0 comments


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

mdsumner avatar Apr 27 '22 18:04 mdsumner