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)