tas_2m_dem.md

isolines_terra_sf <- function(x, levels) {
  if (missing(levels)) {
    levels <- pretty(unlist(minmax(x[[1]])))
    
  }
  ## OMG: note the [[1]] and wide = TRUE which is also weird but different to raster ...
  b <- isoband::isolines(xFromCol(x), yFromRow(x), as.matrix(x[[1]], wide = TRUE), levels = levels)
  
  sf::st_sf(level = levels, geometry  = sf::st_sfc(isoband::iso_to_sfg(b), crs = crs(x)))
}


ex <- c(143, 149, -44, -39)
grat <- as.polygons(rast(ext(ex), res = 1), crs = "OGC:CRS84")
srcdem <- "/vsicurl/https://s3.us-west-2.amazonaws.com/us-west-2.opendata.source.coop/alexgleith/tasmania-dem-2m/Tasmania_Statewide_2m_DEM_14-08-2021.tif"
srccst <- "/vsizip//vsicurl/https://github.com/wmgeolab/geoBoundaries/raw/main/releaseData/CGAZ/geoBoundariesCGAZ_ADM0.zip"  
srcsql <- "SELECT shapeGroup FROM geoBoundariesCGAZ_ADM0 WHERE shapeGroup IN ('AUS')"

library(terra)
tgt0 <- rast(ext(ex), res = 0.002, crs = "OGC:CRS84")

## maybe project to something local (ok to project an empty raster to a string)
tgt0 <- project(tgt0, "EPSG:32755")
## project in case we did with the raster sources
cst <- project(vect(srccst, query = srcsql), tgt0)

dem <- trim(project(rast(srcdem), tgt0, by_util = TRUE))
bathy <- project(rast(sds::gebco()), tgt0, by_util = TRUE)
bathy[bathy > 0] <- NA

bathycont <- isolines_terra_sf(bathy, c(-50))#, -100, -150, -200, -250,-500, -1000, -2000 ))

plot(dem,  col = grey.colors(24, 0, 1), axes = F)
plot(cst, add = TRUE, border = "green")
plot(bathycont, add = TRUE, col = "firebrick")
plot(project(grat, tgt0), add = T)