The end of the world ... how to deal with it?

[mtennekes]

One of the very few features in v3 that is not (yet) in v4:

tm_shape(World, projection = "+proj=eck4") + tm_polygons() + tm_style("natural")

This "natural" style has a tmap option earth.boundary set to TRUE.

The v3 implementation is beyond ugly, so I'm looking for a more elegant solution. The main problem was/is that the 'boundary' of the earth is not simply the longlat bounding of -180 - 180 and -90 - 90. A crs with the center at the Pacific or at one of the poles has different earth boundaries.

Does someone know an elegant solution to find the outline of a crs? @Nowosad @edzer @mdsumner

[edzer]

Is this not an issue that could/would be solved at the level of PROJ?

[mdsumner]

gdal_footprint looks pretty good

library(terra)
r <- project(setValues(rast(), 1), "+proj=eck4")
writeRaster(r, tf <- tempfile(fileext = ".tif"), overwrite = TRUE)
system(sprintf("gdal_footprint %s %s", tf, out <- tempfile(fileext = ".fgb")))

plot(vect(out))

I tried running that with the osgeo.gdal api but don't really have the fortitude for that rn.

[mtennekes]

Nice pragmatic approach @mdsumner

In the PROJ issue the function proj_trans_bounds was suggested. Is this function available or could this become available in one of the r-spatial packages (sf, PROJ, ...)?

[mdsumner]

Certainly in PROJ, I'll have a try at some point, but don't have much capacity atm.

[edzer]

See https://github.com/r-spatial/sf/issues/2415

> st_bbox(st_as_stars()) |> st_transform("+proj=eck4")
     xmin      ymin      xmax      ymax 
-16921203  -8460601  16921203   8460601 
...
> st_bbox(vect(out))
     xmin      ymin      xmax      ymax 
-16921203  -8424632  16908719   8460601 

Or the longer path:

> st_bbox(st_as_stars()) |> 
    st_as_sfc() |> 
    st_set_crs(NA) |> 
    st_segmentize(1) |> 
    st_set_crs("OGC:CRS84") |> 
    st_transform("+proj=eck4") -> x
> st_bbox(x)
     xmin      ymin      xmax      ymax 
-16921203  -8460601  16921203   8460601 

where the hoop through NA_crs_ is needed to get the zero distance lines at the poles.

[mdsumner]

I didn't think that was what was being asked (right?), and I think Even misunderstood too. I was surprised to see that proj_trans_bounds returns just the 4 numbers. Confirmed with

library(reticulate)
pyproj <- import("pyproj")

## I have to fix this somewhere back up
pyproj$datadir$set_data_dir("/usr/local/share/proj")
trans <- Transformer$from_crs("EPSG:4326", "+proj=eck4")
trans$transform_bounds(left = -180, bottom = -90, right = 180, top = 90, direction = pyproj$enums$TransformDirection$FORWARD)

[[1]]
[1] -8460601

[[2]]
[1] -8323299

[[3]]
[1] 8460601

[[4]]
[1] 8323299

edit: heck though the numbers are different, unsure ...

It's the same as raster::projectExtent and reproj::reproj_extent, it's not the visible region boundary, just a fleshed out projected extent. fwiw, I use this every day, it's impossible to exist without it.

reproj::reproj_extent(c(-180, 180, -90, 90), "+proj=eck4", source = "EPSG:4326")
[1] -16921203  16921203  -8460601   8460601

[mtennekes]

Yes, that is correct @mdsumner : I'm looking for the visible boundary, expecting a polygon rather than a box.

[edzer]

st_bbox(st_as_stars()) |> 
    st_as_sfc() |> 
    st_set_crs(NA) |> 
    st_segmentize(1) |> 
    st_set_crs("OGC:CRS84") |> 
    st_transform("+proj=eck4") |> 
    plot(axes = TRUE)

[mdsumner]

that only works in this particular case where the normal boundary aligns with the target one, not every projected region is the outline of -180,180,-90,90 and it won't work for global laea at all (as just one example)

[mdsumner]

also @mtennekes I can't see how to distinguish the edge of the world from NA, in your above example there are missing values literally outside the boundary which is convenient for the footprint approach, but here for example the NAs are just arbitrarily applied

library(terra)
r <- rast("/vsicurl/https://noaadata.apps.nsidc.org/NOAA/G02135/south/daily/geotiff/2024/08_Aug/S_20240823_concentration_v3.0.tif")
pts <- reproj::reproj_xy(xyFromCell(r, 1:ncell(r)), "EPSG:4326", source = crs(r))

r[pts[,2] > -50] <- NA
writeRaster(r, tf <- tempfile(fileext = ".tif"), overwrite = TRUE)
system(sprintf("gdal_footprint %s %s", tf, out <- tempfile(fileext = ".fgb")))

plot(vect(out))

and, there are many projections that will simply repeat (it goes out to some factor of pi I noticed, and it's not symmetric in each direction but I'm not sure). So I'm not sure the question is entirely well defined.

this for example, it's nice that the warper now does this, but I think you're asking for a boundary that doesn't always exist, although I think you could reasonably expect something. 🙏

library(terra)
dsn <- "/vsicurl/https://gebco2023.s3.valeria.science/gebco_2023_land_cog.tif"
r <- project(rast(dsn), rast(ext(c(-1, 1, -1, 1) * 5e7), ncols = 1024, nrows = 1024, crs = "+proj=omerc +lonc=147 +alpha=15"), by_util = TRUE)
plot(r)

[mtennekes]

st_bbox(st_as_stars()) |> 
    st_as_sfc() |> 
    st_set_crs(NA) |> 
    st_segmentize(1) |> 
    st_set_crs("OGC:CRS84") |> 
    st_transform("+proj=eck4") |> 
    plot(axes = TRUE)

Thx @edzer Something along these lines would be sufficient for compatibility with tmap3. However, a test with EPSG 3035 wasn't successful yet:

EU = st_transform(World, crs = 3035)

poly = st_bbox(st_as_stars()) |> 
	st_as_sfc() |> 
	st_set_crs(NA) |> 
	st_segmentize(1) |> 
	st_set_crs("OGC:CRS84") |> 
	st_transform(3035) |> 
	st_as_sf()

tm_shape(EU) + 
	tm_polygons() +
	tm_shape(poly) +
	tm_borders(col = "purple", lwd = 3)

Yes @mdsumner I'm fully aware that not every crs has such a visible boundary (not even WGS84/Mercator)

[DOSull]

Yes @mdsumner I'm fully aware that not every crs has such a visible boundary (not even WGS84/Mercator)

Especially not Mercator!

A tmap function that resolved this problem for even a fraction of common world projections would be a major achievement as it's very hard problem. A solution requires inferring for any given projection what it's boundary conditions look like. There simply isn't going to be a one-size fits all answer for that problem, at any rate not a simple one at the level I assume you are hoping for.

Finite cyclindrical rectangular projections or even pseudo-cylindrical projections, if you know their central meridian, would allow a solution similar to @edzer's with the great circle appropriately modified, but in many cases the underlying projection will be making such a mess of the world map that colouring in the 'world footprint' appropriately is the least of your problems!

World$geometry |> st_transform("+proj=moll +lon_0=120") |> plot()

That's a simple central meridian shift. Once you start working with oblique forms and shifting projection centres things get even messier. It's a problem that d3 does a much better job of than proj, but that implies a major overhaul of proj and even how we approach geospatial data!

If you do figure out an approach, sing it from the rooftops!

[edzer]

An attempt to cut out the visible part for orthographic projections is found here, but I'd admit it's pretty ugly.

[mtennekes]

I've integrated Edzers solution

st_bbox(st_as_stars()) |> 
    st_as_sfc() |> 
    st_set_crs(NA) |> 
    st_segmentize(1) |> 
    st_set_crs("OGC:CRS84") |> 
    st_transform("+proj=eck4") |> 
    plot(axes = TRUE)

into tmap4, which works well for many standard pseudo-cyclindrical projections. That makes it backwards compatible with tmap3, which is in this stadium sufficient. Later, we can see if other projections families (e.g. orthogonal) can be supported as well.

world_projs = 
  c("wintri", 
     "robin", 
    paste0("eck", 1:6), 
    "moll", 
    paste0("wag", 1:7), "eqearth")

tms = lapply(world_projs, function(p) {
	tm_shape(World, crs = paste0("+proj=", p)) + tm_polygons() + tm_style("natural") + tm_title(p)
})
do.call(tmap_arrange, c(tms, list(ncol= 4)))