graticule_map.md

library(graticule)
lons <- seq(-80, -40, by = 5)
lats <- seq(77, 84.5, by = .5)
latpts <- cbind(min(lons), lats)
lonpts <- cbind(lons, min(lats))

prj <- "+proj=laea +lat_0=90 +lon_0=-60"
xy_lon <- sf::sf_project(pts = lonpts, to = prj, from = "OGC:CRS84")
xy_lat <- sf::sf_project(pts = latpts, to = prj, from = "OGC:CRS84")
g <- graticule(lons,lats, proj = prj)
plot(g)

lonlab <- function(x) paste(x, ifelse(x >0, "E", "W"))
latlab <- function(x) paste(x, ifelse(x >0, "N", "S"))

text(xy_lon, lab = lonlab(lonpts[,1]), pos = 1, offset = 1)
text(xy_lat, lab = latlab(latpts[,2]), pos = 1, offset = -4)

get the ice data

  library(graticule)
#> Loading required package: sp
lons <- seq(-80, -40, by = 5)
lats <- seq(77, 84.5, by = 1)
latpts <- cbind(min(lons), lats)
lonpts <- cbind(lons, min(lats))

#prj <- "EPSG:3413"
prj <- "+proj=stere +lon_0=-60 +lat_0=90"
xy_lon <- sf::sf_project(pts = lonpts, to = prj, from = "OGC:CRS84")
xy_lat <- sf::sf_project(pts = latpts, to = prj, from = "OGC:CRS84")
pts <- rbind(xy_lon, xy_lat)
ex <- c(range(pts[,1]), range(pts[,2]))
g <- graticule(lons,lats, proj = prj, tiles = TRUE)
sf::write_sf(sf::st_as_sf(g), gpkg <- tempfile(fileext = ".gpkg"))
#elevsrc <- "/vsicurl/https://public.services.aad.gov.au/datasets/science/GEBCO_2021_GEOTIFF/GEBCO_2021.tif"
icesrc <- "/vsicurl/https://seaice.uni-bremen.de/data/amsr2/asi_daygrid_swath/n3125/2019/aug/Arctic3125/asi-AMSR2-n3125-20190830-v5.4.tif"
ice <- whatarelief::elevation(source = icesrc, extent = ex, projection = prj, dimension = dev.size("px"), 
                              options = c("-cutline", gpkg), resample = "bilinear")

ice[!ice > 0] <- NA
ice[ice > 100] <- NA

tf <- tempfile(fileext = ".tif")
im <- whatarelief::elevation(extent = ex, projection = prj, dimension = dev.size("px"), 
                              options = c("-cutline", gpkg))
#> [1] "/vsicurl/https://public.services.aad.gov.au/datasets/science/GEBCO_2021_GEOTIFF/GEBCO_2021.tif"

par(xpd = NA)
ximage::ximage(ice, extent = ex, asp = 1, axes = FALSE, col = hcl.colors(24))
ximage::xcontour(im, extent = ex, add = TRUE, levels = 0)
plot(g, add = T)

lonlab <- function(x) paste(x, ifelse(x >0, "E", "W"))
latlab <- function(x) paste(x, ifelse(x >0, "N", "S"))

text(xy_lon, lab = lonlab(lonpts[,1]), pos = 1, offset = 1, cex = 0.7)
text(xy_lat, lab = latlab(latpts[,2]),  pos = 4, offset = -2, cex = 0.7)

^{Created on 2023-01-12 by the reprex package (v2.0.1)}

Excellent! Let me process this and I will get back :)

What is the primary use of whatarelief? I really need to follow more the work you are doing :)

all good, this prompted me to push a few necessary changes forward 👌 we'll do this as examples with various options of how 🙏

💯

whatarelief is a helper experiment around a pure- GDAL reader, in between pure-read and visualisation- or format-requirements, it's about having componentry tools not everything bundled into one do-all monster package

(sf tries to dominate everything, but really fails miserably when you step even remotely off the path you're meant to stick to, I can bend it where I want to go but only after exploring the parts properly 😇)

I hear you :) I will certainly have a look at whatarelief soon and integrate it to my workflow.