marginal polygons missing when plotted with orthographic projection

[AMBarbosa]

I'm not sure if this is intended or unavoidable behaviour, but when plotting a vector map in orthographic projection, polygons that are partly outside the mapping region are entirely left out. This can be quite conspicuous for large polygons, e.g. mainland North America, Russia and Antarctica in the following map:

wrld <- geodata::world(path = tempdir())

earth_centr <- vect(cbind(0, 0))
earth <- buffer(earth_centr, width = 6371000)  # Earth radius

ortho <- "+proj=ortho +lat_0=0 +lon_0=-20"

plot(earth, col = "lightblue", border = NA, axes = FALSE)
plot(project(wrld, ortho), col = "tan", add = TRUE)

Warning messages: 1: Point outside of projection domain (GDAL error 1)

Would there be a way to crop the vector map to the mapped / visible region before plotting, to avoid throwing away entire polygons? Cheers!

[mdsumner]

I tried this and it seems ok (basically use your disc at a slightly narrower radius to crop the source in the longlat crs). It's often the case that you can find a solution like this for a given dataset with a bit of tweaking. General solutions remain elusive (though see the work in D3 by Mike Bostock: https://bost.ocks.org/mike/example/ - sadly the graphics seem broken atm but it's a really excellent tale).

## make another a bit smaller
earth1 <- buffer(earth_centr, width = 6350000)  

set.crs(earth1, ortho)
## then crop the source in original crs
xx <- project(earth1, crs(wrld))
plot(earth, col = "lightblue", border = NA, axes = FALSE)

plot(project(crop(wrld, xx), ortho), col = "tan", add = TRUE)

[rhijmans]

I added argument "partial=FALSE" to project. Set it to TRUE to keep geometries that can only partially be represented in the target CRS. That helps:

library(terra)
wrld <- geodata::world(path =".")
earth_centr <- vect(cbind(0, 0))
earth <- buffer(earth_centr, width = 6371000)  # Earth radius
ortho <- "+proj=ortho +lat_0=0 +lon_0=-20"

p <- project(wrld, ortho, partial=TRUE)
plot(earth, col = "lightblue", border = NA, axes = FALSE)
plot(p, col = "tan", add = TRUE)

And this also works when centering the projection near the dateline. But edges remain a problem because there are insufficient nodes there. That can be fixed like this.

p <- project(p, wrld)
w <- crop(wrld, ext(p))
w <- densify(w, 100000)
pw <- project(w, ortho)

plot(earth, col = "lightblue", border = NA, axes = FALSE)
plot(pw, col = "tan", add = TRUE)

But that does not always work for Antarctica (densify has no effect because at the South Pole distances are zero).

Projection to "ortho" probably needs a specialized treatment.

[mdsumner]

crop it to -89.9 or similar is a common trick, a bit like the disc-project-crop (and also works for the -180/180 ambiguity when zingers occur)

but, can't densify work independent of crs units? I got frustrated by this in sf, and in my ideal world there would be functions to work with numbers as-is, and functions to work with unit/geo-awareness - but I had zero impact on where sf went with that

[rhijmans]

@mdsumner: I forgot, but you can set use argument "flat=TRUE" in densify to work with the raw numbers.

Here is another attempt, using the current terra-dev (for "crop" that wraps lon/lat extents around the date-line). This seems to work OK (but only for lat_0=0).

library(terra)
wrld <- geodata::world(path =".")
earth_centr <- vect(cbind(0, 0))
earth <- buffer(earth_centr, width = 6371000)

ortho <- function(x, lon) {
    nlon <- lon - 90
    xlon <- lon + 90
    v <- crop(x, c(nlon, xlon, -90, 90), TRUE)
    prj <- paste0("+proj=ortho +lat_0=0 +lon_0=", lon)
    v <- densify(v, 1, flat=TRUE)
    project(v, prj, partial=TRUE)
}

plortho <- function(lon) {
    p <- ortho(wrld, lon)
    plot(earth, col = "lightblue", border = NA, axes = FALSE)
    plot(p, col = "tan", add = TRUE)
}

plortho(-80)
plortho(20)
plortho(100)

[mdsumner]

awesome, thank you 🙏

[mdsumner]

to my mind the logic should be a function for ellipsoid aware logic, and a function for flat logic, no matter what projection the source is in. it bugs me that it depends on the input, flat = F will have no impact for projected data. They are different functions and I don't get why it pivots on crs type. But, terra is way more sensible than alternatives and I'm v grateful for that

[rhijmans]

@mdsumner I think that when computing the area of polygons it is reasonable to assume that one wants to output in m2 or such, not in "square degrees". But if you really want that, you can assign a planar CRS to your SpatVector. That is annoying and therefore in a function like densify, where it seems more likely that you might want to treat the lon/lat coordinates as planar there is a little helper in "flat=TRUE". The alternative interface that you suggest, separate functions, would be cleaner, but then you would need to decide which function to use to, for example, compute area of polygons, by first inspecting the CRS. I think that would lead to a lot of mistakes and clumsy code; so I do not prefer it.

[mdsumner]

I don't think it is reasonable, if your data is in a projection you get "square metres" but big deal, it depends where in what crs if that is close to the real area. I want the "area" that the data has. I I want the area on the ground, I use an appropriate method. If it's in Mercator or Stereographic or so many others ... the relationship of the data area to the real world depends where in the domain the measurement is and the properties of the crs. It's a specious trap for non experts I don't understand the appeal. Id like a "on the ground area" calculator to work by always using ellipsoidal methods ( unproject, calculate for the Earth, etc). The area area can just be what you see on the screen, this current situation is much more confused, literally specious.

Anyway, as I said, this appeal had zero impact elsewhere... but my ongoing advice will always be don't trust these tools, you have to understand the details.

[mdsumner]

oh wait, I'm sorry - expanse does exactly this by default - wow, that's good

I'm not across what terra is doing enough, will explore more

[mdsumner]

love it, so glad to see this made clear:

For vector data, the best way to compute area is to use the longitude/latitude CRS. This is contrary to (erroneous) popular belief that suggest that you should use a planar coordinate reference system. This is done automatically, if transform=TRUE.