hypertidy / gdalio

consider mdsumner/whatarelief instead WIP (helpers for gdal warper, and wrappers for various packages)
Other
27 stars 1 forks source link

stac examples #19

Open mdsumner opened 2 years ago

mdsumner commented 2 years ago

just a place to store some examples while I figure out this stac thing, shared here for @cboettig

## goal, convert this python notebook to R
#https://planetarycomputer.microsoft.com/dataset/naip#Example-Notebook
## discussion ongoing of Michael Sumner and Carl Boettiger

ext <- c(-111.9839859008789, -111.90502166748045, 40.5389819819361, 40.57015381856105)

range_old = "2010-01-01/2013-01-01"
range_new = "2018-01-01/2020-01-01"
library(rstac)
catalog <- stac("https://planetarycomputer.microsoft.com/api/stac/v1")

search_old <- catalog %>%  stac_search(
    collections="naip", bbox=ext[c(1, 3, 2, 4)], datetime=range_old
)
search_new = catalog %>% stac_search(
    collections="naip", bbox  = ext[c(1, 3, 2, 4)], datetime=range_new
)

item_old <- search_old %>% get_request()
item_new <- search_new %>% get_request()
url_old <- sprintf("/vsicurl/%s", item_old$features[[1]]$assets$image$href)
url_new <- sprintf("/vsicurl/%s", item_new$features[[1]]$assets$image$href)

## these can be slow, sometimes it takes a while
## we just getting 'raster info', like gdalinfo output
## (it makes our reprex very slow ... but sometimes it's fast ...)
ri_old <- vapour::vapour_raster_info(url_old) 
ri_new <- vapour::vapour_raster_info(url_new)

## we can use a longlat grid like they do in the notebook (guessing at dim here ...)
#gdalio_set_default_grid(list(extent = ext, dimension = c(512, 512), projection = "OGC:CRS84"))
## or, we can get the actual grid spec and use it (but use a lower res overview - there are 5 levels in total, 
## $dimXY and then 4 more in $overviews)

## new and old are entirely different, but the warper doesn't care we'll read both to match the second
## lowest overview of old
library(gdalio)
gdalio_set_default_grid(list(extent = ri_old$extent, 
                             dimension = tail(ri_old$overviews, 4)[1:2], 
                             projection = ri_old$projection))

## load code for terra etc formats
source(system.file("raster_format/raster_format.codeR", package = "gdalio", mustWork = TRUE))
## we may get obscure terra error/warnings but that's 'normal'
image_old <- gdalio_terra(url_old, bands = 1:3)
image_new <- gdalio_terra(url_new, bands = 1:3)
op <- par(mfrow = c(1, 2))
library(terra)
#> terra version 1.4.17
#> 
#> Attaching package: 'terra'
#> The following object is masked from 'package:gdalio':
#> 
#>     vrt
plotRGB(image_old)
plotRGB(image_new)

par(op)

Created on 2021-12-04 by the reprex package (v2.0.1)

mdsumner commented 2 years ago

@cboettig this one you shared in slack, it doesn't work for me ??

library(rstac)
  s_obj <- stac("https://planetarycomputer.microsoft.com/api/stac/v1/")

it_obj <- s_obj %>% 
  stac_search(collections = "landsat-8-c2-l2",
              bbox = c(-47.02148, -17.35063, -42.53906, -12.98314)) %>%
  get_request() %>%
  items_sign(sign_fn = sign_planetary_computer())
#> Warning: Items matched not provided.
#> Error: No collection found with id 'landsat-8-c2-l2'

Created on 2021-12-04 by the reprex package (v2.0.1)

cboettig commented 2 years ago

weird, I think I recall some network instability sometimes not getting me results. Still working ok for me now:

library(rstac)
s_obj <- stac("https://planetarycomputer.microsoft.com/api/stac/v1/")

it_obj <- s_obj %>% 
  stac_search(collections = "landsat-8-c2-l2",
              bbox = c(-47.02148, -17.35063, -42.53906, -12.98314)) %>%
  get_request() %>%
  items_sign(sign_fn = sign_planetary_computer())
#> Warning: Items matched not provided.
it_obj
#> ###STACItemCollection
#> - features (250 item(s)):
#>   - LC08_L2SP_219072_20211126_02_T1
#>   - LC08_L2SP_219071_20211126_02_T1
#>   - LC08_L2SP_219070_20211126_02_T1
#>   - LC08_L2SP_219069_20211126_02_T1
#>   - LC08_L2SP_221071_20211124_02_T1
#>   - LC08_L2SP_221070_20211124_02_T1
#>   - LC08_L2SP_221069_20211124_02_T1
#>   - LC08_L2SP_218072_20211119_02_T1
#>   - LC08_L2SP_218071_20211119_02_T1
#>   - LC08_L2SP_218070_20211119_02_T1
#>   - ... with 240 more feature(s).
#> - assets: 
#> ANG, SR_B1, SR_B2, SR_B3, SR_B4, SR_B5, SR_B6, SR_B7, ST_QA, ST_B10, MTL.txt, MTL.xml, ST_DRAD, ST_EMIS, ST_EMSD, ST_TRAD, ST_URAD, MTL.json, QA_PIXEL, ST_ATRAN, ST_CDIST, QA_RADSAT, thumbnail, SR_QA_AEROSOL, reduced_resolution_browse, tilejson, rendered_preview
#> - other field(s): type, features, links

Created on 2021-12-04 by the reprex package (v2.0.1)

sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C              LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] arrow_6.0.0.20211112 raster_3.5-2         sp_1.4-6             rstac_0.9.1-5        dplyr_1.0.7          testthat_3.1.0      
[7] terra_1.4-22        

loaded via a namespace (and not attached):
 [1] styler_1.6.2      progress_1.2.2    tidyselect_1.1.1  xfun_0.28         purrr_0.3.4       lattice_0.20-44   vctrs_0.3.8      
 [8] generics_0.1.1    htmltools_0.5.2   yaml_2.2.1        utf8_1.2.2        rlang_0.4.12      R.oo_1.24.0       pillar_1.6.4     
[15] glue_1.5.1        withr_2.4.3       DBI_1.1.1         R.utils_2.11.0    bit64_4.0.5       R.cache_0.15.0    lifecycle_1.0.1  
[22] R.methodsS3_1.8.1 zip_2.2.0         bench_1.1.1       codetools_0.2-18  evaluate_0.14     knitr_1.36        callr_3.7.0      
[29] tzdb_0.2.0        fastmap_1.1.0     ps_1.6.0          curl_4.3.2        fansi_0.5.0       highr_0.9         Rcpp_1.0.7       
[36] readr_2.1.0       backports_1.3.0   desc_1.4.0        pkgload_1.2.3     vroom_1.5.6       jsonlite_1.7.2    fs_1.5.0         
[43] bit_4.0.4         hms_1.1.1         digest_0.6.29     processx_3.5.2    duckdb_0.3.1-1    grid_4.1.0        rprojroot_2.0.2  
[50] cli_3.1.0         tools_4.1.0       magrittr_2.0.1    tibble_3.1.6      crayon_1.4.2      neonstore_0.4.4   pkgconfig_2.0.3  
[57] ellipsis_0.3.2    prettyunits_1.1.1 reprex_2.0.1      rstudioapi_0.13   assertthat_0.2.1  rmarkdown_2.11    httr_1.4.2       
[64] R6_2.5.1          compiler_4.1.0