Wanted: high(er) res bathymetry option

[Maschette]

Similar to rnatural earth where you can specify the scale/size/resolution of the polygons for countries it would be great to have an option for higher res bathymetry. The current one is great for the full SOmap() but it is quite blocky for small scale areas using SOmap_auto.

One option would be have an extra package that has a few resolutions of bathy and uses those when told to.

[raymondben]

I'm reluctant, TBH. The data size in this package is already excessive, so bundling more is probably not advisable. A separate package is an option, as you say - but we already have packages that facilitate access to bathy and similar gridded data. And users can pass their own bathy raster to SOmap_auto already ... so ... not a priority by my thinking.

If the issue is primarily about the blocky nature of the visual presentation, rather than a genuine need for higher res data (e.g. as you might want for analysis purposes rather than illustration purposes) then perhaps an alternative might be to cook up an example for users that shows how to convert the bathy raster into contours, which will appear suitably smooth (and won't be too horrendously inaccurate, at least at moderate spatial scales)?

[Maschette]

This just popped into my head after talking with someone else who wanted a higher res bathy in the background. @mdsumner do you still have the code you used to make the initial 16km grid?

[mdsumner]

of course:

https://github.com/AustralianAntarcticDivision/SOmap/blob/master/data-raw/Bathy.R

That used GEBCO 2014 and now we have 2019, with those (huge) files on hand it's pretty easy to create a projected raster at whatever resolution is required - I might use stars or terra now but it won't make a massive difference

[mdsumner]

Walkthrough the code (because it's weird)

## the map projection, this is the default for SOmap (technically we should drop ellps etc)
prj <- "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m"   # +no_defs +ellps=WGS84 +towgs84=0,0,0

## read the whole thing into memory (because I can)
src <- readAll(raadtools::readtopo("gebco_14"))  ##, xylim = extent(-180, 180, -90, 0))

proj_longlat <- "+proj=longlat +datum=WGS84"

## specify the target grid, this is a shortcut way to get a projection extent from 
## the southern ocean 90S-5S
r <- raster(projectExtent(raster(extent(-180, 180, -90, -5), crs  = proj_longlat), prj))
## but I want clean extents and clean pixel sizes at 16000 on whole numbers so
r <- raster(spex::buffer_extent(r, 16000), crs = prj)
res(r) <- 16000
## now do the hard work (might be impractical with GEBCO2019, but there's ways)
Bathy <- projectRaster(src, r)
## set the smallest size possible but still keep sense of the values (because data in pkg)
## this is 16 bit integer which is enough, signed 
dataType(Bathy) <- "INT2S"
## crush the values into integer and create the final result for saving to file
Bathy <- setValues(Bathy, as.integer(values(Bathy)))

[mdsumner]

There's tools like piggyback to take a data set like this and store it in our repo here, then a function could download it without sullying the perfect tiny size of SOmap in its pure form. ;)

[Maschette]

[mdsumner]

yeah it's not trivial, gebco_19 is massive overkill I shouldn't have mentioned it - shall we go for 8km and 4km?

[mdsumner]

i.e.

prj <- "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m"   # +no_defs +ellps=WGS84 +towgs84=0,0,0
library(raadtools)
src <- readAll(raadtools::readtopo("gebco_14", xylim = extent(-180, 180, -90, 0)))

proj_longlat <- "+proj=longlat +datum=WGS84"

## specify the target grid
r <- raster(SOmap::Bathy)
grain <- 8000
res(r) <- grain
## now do the hard work
Bathy <- projectRaster(src, r)

But, you might also encourage folks to get very local topography using ceramic

i.e.

library(ceramic)
library(raster)
elev <- cc_location(extent(100, 120, -70, -50), type = "elevation-tiles-prod")

That stuff comes out in Mercator, so could be projected to whatever target's in use - it picks a zoom based on the region asked for, a kind of rough number of tiles it has to download (they get cached).

[Maschette]

I was thinking 8km, 4km and possibly 2km dependent on size in a separate package. Ideally I would like to use gebco19 as the base that we grid off as in some regions the depths are quite different but it would probably be a pain in the ass to make.

I think we could definitely write a guide for using ceramic as the bathymetry as well.

[mdsumner]

ok I've done something in https://github.com/AustralianAntarcticDivision/SOmap/pull/93

but it involves a lot of downloads and stars and raster and file churn ...

a proper fix in the works from the GDAL warp stuff I've been working on

[mdsumner]

ok done in https://github.com/AustralianAntarcticDivision/SOmap/pull/96

try it out @Maschette

[mdsumner]

some examples, will work anywhere

library(SOmap)
#> Loading required package: raster
#> Loading required package: sp
SOmap_auto(c(143, 149), c(-46, -41))
#> [1] 512 512

SOmap_auto(c(80, 90), c(-64, -58))
#> [1] 512 512
#> Warning in image.default(x = x, y = y, z = value, useRaster = useRaster, :
#> unsorted 'breaks' will be sorted before use

SOmap_auto(c(153, 164), c(-58, -54), dimXY = dev.size("px"), input_lines = FALSE)
#> [1] 480 672
#> Warning in image.default(x = x, y = y, z = value, useRaster = useRaster, :
#> unsorted 'breaks' will be sorted before use

^{Created on 2021-10-13 by the reprex package (v2.0.1)}

[Maschette]

So how does this work if you are offline?

[mdsumner]

it'll fall back to the old, but you could insert <path to local topo file> in place of the online source in .get_raad_gebco()

[Maschette]

Awesome, are the warnings from image.default needed?

[mdsumner]

what warnings? can you please report details and have examples?

I'll put a tentative patch into branch get-bathy-gebco-raad to allow specifying a local topo file

[Maschette]

In your example above there are warnings after the latter 2:

SOmap_auto(c(153, 164), c(-58, -54), dimXY = dev.size("px"), input_lines = FALSE)
#> [1] 480 672
#> Warning in image.default(x = x, y = y, z = value, useRaster = useRaster, :
#> unsorted 'breaks' will be sorted before use

[mdsumner]

oh good ... in my reprex ... I have no idea please stick to the topic and follow up issues independently, I don't see that having any impact on the current issue

[mdsumner]

this is delivered, more to come in doc etc

[Maschette]

oh good ... in my reprex ... I have no idea please stick to the topic and follow up issues independently, I don't see that having any impact on the current issue

Thanks Mike, I was not sure if it was a result of this change or not. After some checking it is not so I will try and track it down.

[Maschette]

Also, I forgot to mention, fucking great work on this. It looks awesome.