cutline examples

[mdsumner]

worked on a local file

d <- vapour:::gdal_raster_data(whatarelief::gebco22(), target_dim = dev.size("px"),  options = c("-cutline", system.file("gpkg/nc.gpkg", package = "sf")))

but that code crashes R

[mdsumner]

fine at the command line

export cog=/vsicurl/https://gebco2022.s3.valeria.science/gebco_2022_complete_cog.tif
gdalwarp $cog out.tif -cutline autotest/ogr/data/testfgb/wfs/get_feature.fgb -cl unknown -ts 100 100 -crop_to_cutline

[mdsumner]

well it does work, so we just might need a heap more testing - this was in 0.9.3.9005

R 4.2.3; x86_64-pc-linux-gnu; 2023-03-27 20:52:27 UTC; unix

crash was in slightly later main

curl::curl_download("https://github.com/r-spatial/sf/raw/main/inst/gpkg/nc.gpkg", "nc.gpkg")

cog <- "/vsicurl/https://gebco2022.s3.valeria.science/gebco_2022_complete_cog.tif"

vec <- "nc.gpkg"

v <- vapour:::gdal_raster_data(cog, options = c("-cutline", vec, "-crop_to_cutline", "-cl", "nc.gpkg"), target_dim = c(100, 100))

[mdsumner]

hell yeah, here's using from purely online sources

vurl <- "/vsicurl/https://github.com/r-spatial/sf/raw/main/inst/gpkg/nc.gpkg"
cog <- "/vsicurl/https://gebco2022.s3.valeria.science/gebco_2022_complete_cog.tif"
layers <- vapour::vapour_layer_names(vurl)
v <- vapour:::gdal_raster_data(cog, options = c("-cutline", vurl, 
                                                "-crop_to_cutline", 
                                                "-cl", layers[1L]), target_dim = c(100, 100))

f <- vapour:::gdal_raster_dsn(cog, options = c("-cutline", vurl, 
                                               "-crop_to_cutline", 
                                               "-cl", layers[1L]), target_dim = c(100, 100))

str(v)
List of 1
 $ : num [1:10000] 0 0 0 0 0 0 0 0 0 0 ...
 - attr(*, "dimension")= num [1:2] 100 100
 - attr(*, "extent")= num [1:4] -84.3 -75.5 33.9 36.6
 - attr(*, "projection")= chr "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AU"| __truncated__
> str(f)
List of 1
 $ : chr "/perm_storage/home/data/r_tmp//RtmpYJROR9/file3dbe227bd6d9f.tif"
 - attr(*, "dimension")= num [1:2] 100 100
 - attr(*, "extent")= num [1:4] -84.3 -75.5 33.9 36.6
 - attr(*, "projection")= chr "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AU"| __truncated__

[mdsumner]

more examples

it's neat having a SQL clause to pick features

dsn <- "/vsizip//vsicurl/https://github.com/wmgeolab/geoBoundaries/raw/main/releaseData/CGAZ/geoBoundariesCGAZ_ADM0.zip"
#where <- c("Australia", "New Zealand", "Antarctica")[1]
where <- "Australia"
codes <- countrycode::countrycode(where, "country.name", "iso3c")

layer <- vapour::vapour_layer_names(dsn)[1L]

csql <- sprintf("SELECT shapeGroup FROM %s WHERE shapeGroup IN (%s)", layer,  paste(paste0("'", codes, "'"), collapse = ","))
g <- vapour::vapour_read_geometry(dsn, sql = csql)
geb <- "/vsicurl/https://public.services.aad.gov.au/datasets/science/GEBCO_2021_GEOTIFF/GEBCO_2021.tif"
d <- vapour::gdal_raster_data(geb, target_dim = c(100, 100), options = c("-cutline", dsn, "-csql", csql, "-crop_to_cutline"))
whatarelief:::imfun(d, coastline = F)

where <- "Antarctica"
codes <- countrycode::countrycode(where, "country.name", "iso3c")

csql <- sprintf("SELECT shapeGroup FROM %s WHERE shapeGroup IN (%s)", layer,  paste(paste0("'", codes, "'"), collapse = ","))
src <- sprintf("vrt://%s?scale=0,255", whatarelief:::.imagery_sources["wms_virtualearth"])
d <- vapour::gdal_raster_data(geb, target_dim = c(512, 512)* 2, target_crs = "EPSG:3031",
                              options = c("-cutline", dsn, "-csql", csql, "-crop_to_cutline"))
whatarelief:::imfun(d, coastline = F)