Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Last active August 21, 2021 11:55
Show Gist options
  • Save mdsumner/fe83a9b1e6e0d49bf449ef40df690dd6 to your computer and use it in GitHub Desktop.
Save mdsumner/fe83a9b1e6e0d49bf449ef40df690dd6 to your computer and use it in GitHub Desktop.
stars example (with raster help) to build the read_stars RasterIO list from an extent
  f <- "ftp.nass.usda.gov/download/res/2019_30m_cdls.img"
  library(raster)
#> Loading required package: sp
crop_io <- function(x, ext, resample = "nearest_neighbour") {
  rx <- raster::raster(x)
  if (inherits(ext, "bbox")) {
    ext <- raster::extent(ext[c("xmin", "xmax", "ymin", "ymax")])
  } else {
    ext <- raster::extent(ext)
  }
  col1 <- raster::colFromX(rx, ext@xmin)
  col2 <- raster::colFromX(rx, ext@xmax)
  row1 <- raster::rowFromY(rx, ext@ymax)
  row2 <- raster::rowFromY(rx, ext@ymin)
  if (is.na(row1)) row1 <- 1
  nx <- col2 - col1
  ny <- row2 - row1
  if (is.na(row2)) row2 <- rx@nrows
  
                   list(nXOff = col1,
                        nYOff = row1,
                        nXSize = nx,
                        nYSize = ny,
                        nBufXSize = nx,
                        nBufYSize = ny,
                        resample = resample)                   
}
r0 <- raster::raster(f)
index <- raster::raster(raster::extent(r0), 
                        res = raster::res(r0)*512)

st0 <- stars::read_stars(f, proxy = TRUE)
data("wrld_simpl", package = "maptools")
sfx <- sf::st_transform(sf::st_as_sf(subset(wrld_simpl, NAME == "United States"))
                        , sf::st_crs(st0))

plot(extent(index))
plot(sfx, add = TRUE)
#> Warning in plot.sf(sfx, add = TRUE): ignoring all but the first attribute

cell <- 25  ## first index tile with non-NA
ext <- raster::extentFromCells(index, cell)
io <- crop_io(f, ext)
st <- stars::read_stars(f, RasterIO = io)
plot(st)

## find interactively
#cell <- raster::click(setValues(index, 0), cell = TRUE, n = 1, show   = FALSE)[, "cell"]
cell <- 26054
ext <- raster::extentFromCells(index, cell)
io <- crop_io(f, ext)
st <- stars::read_stars(f, RasterIO = io)
plot(st)

Created on 2020-05-01 by the reprex package (v0.3.0)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment