mapPlot() handles Antarctica poorly

[dankelley]

Quite a lot of projections seem to work quite well lately, but Antarctica is a problem in several, e.g.

library(oce)
data(coastlineWorld)
par(mar=c(3, 3, 2, 1), mgp=c(2, 0.7, 0))
mapPlot(coastlineWorld, col="gray", axes=FALSE, projection="+proj=ortho")

yields as follows. Perhaps there are problems at other spots, but the real eye sore is Antarctica.

[richardsc]

Antarctica is definitely the worst. Here's an example with mollweide:

library(oce)
data(coastlineWorld)
mapPlot(coastlineWorld, col='gray', projection='+proj=moll +lon_0=180')

[dankelley]

I've started issue code, and 616c.out (at end) is indicating that perhaps part of the problem is difficulty in converting from lon-lat to x-y space. Here's what I get for the ortho projection (on which I'll focus to start):

> data.frame(lon=alon[fake], lat=alat[fake], x=xy$x[fake], y=xy$y[fake])
    lon lat             x        y
1   180 -90 -9.916898e-25 -6370997
2   180 -89           Inf      Inf
3   140 -89           Inf      Inf
4   100 -89           Inf      Inf
5    60 -89  9.629270e+04 -6370027
6    20 -89  3.802896e+04 -6370027
7   -20 -89 -3.802896e+04 -6370027
8   -60 -89 -9.629270e+04 -6370027
9  -100 -89           Inf      Inf
10 -140 -89           Inf      Inf
11 -180 -89           Inf      Inf
12 -180 -90 -2.210189e-25 -6370997

so we see that (a) a lot of points cannot be projected and (b) those that do end up at x=0, not scattered along the edge of the globe's outside arc. I'm thinking a solution might be in the deep-down proj.4 code.

[mdsumner]

It needs segmentizing (adding redundant vertices) along the degenerate parallel at -90. Personally very interested in the variety of available solutions to this, the ultimate of which ( in the planar case) is D3 with Bostock's Flawed Example.

The segmentizing algo needs knowledge of the curvature within the target projection, which is why general solutions are rare. But also different data sets have different specifics to worry about.

[mdsumner]

Oops I see you already segmentizing, which suggests it needs to be offset from -90 by a tiny amount to stay in bounds for some projections. Keen to share what you've learned here

[dankelley]

Thanks @mdsumner, it's great to have more eyes on this!

This particular issue is one of several related problems that are likely to have solutions that will be projection-dependent. (Longitudes +90 and -90 are only at the "edge of the earth" if the view is from above the Greenwich meridian ... which is not the view oceanographers tend to want in looking at the Pacific! Also, latitude can be a problem too, if the viewer's "eye" is not above the equator.)

It's quite hard to compute where the "edge of the earth" is on the projected coordinates (i.e. "page" coordinates) because there are a lot of projections that have to be considered. One approach I've considered is to fit a curve to the edge of "visible" region. But that curve would have to be mighty complicated for interrupted projections, for example. (Oceanographers sometimes use them, split along continents, so as to highlight the oceans.)

I don't want a solution that is bound closely to individual projections, though, because then the coding (and debugging) effort scales as N, where N is the number of supported projections. I need a solution that can be generalized.

The Antarctic case is special, because of the fact that pole is inside the polygon. But there is also a related issue that we in oce have taken to calling "UHL" (for "ugly horizontal lines"); see e.g. #656 but note that UHL crop up in linegraphs as well as image graphs. I think my solution in oce::coastlineCut() works well for this, although there is a remaining problem of infinitely-thin dribbles along the edge of the earth that have not yet solve in code, although the methodology seems simple -- I will try doing a pass in the projected space, i.e. the "page" space, removing subregions that contain no area; this will work for coastlines (but not other entities) because no sensible coastline contains regions with zero area. (By "region" I mean portions within a polygon ... readers probably understand what I mean but, if not, it will be simpler for me to code and document the solution than to try to explain it in words here.)

Some solutions methods might be easier if all proj4 projections had (working) inverses. Readers might consult #628 and https://stat.ethz.ch/pipermail/r-sig-mac/2016-May/011931.html for a discussion of issues related to inverses, and especially the 'wintri' projection, which had a nonconverging approximation loop, a couple of years ago. (The advantage of "wintri" is an odd one ... other projections that have clean inverses look nearly the same, but there seems to be some sociological benefit, if not technical benefit, in using the projection that's on the wall in most middle-school classrooms.)

For a while, I had within oce my own version of proj4, in which I had corrected some such problems, but I eventually gave up, and simply disallowed some projections in 'oce'. Also, of course, proj4 is a work in progress and it keeps improving, and it seems preferable for people like me to point out problems (and maybe solutions) in proj4, so benefits can be accrued across a wide range of descendent codes.

Again, thanks for the comments.

[mdsumner]

Great, this is really helpful! I want to bring all this out into the open and have better tools for it in R. It's long overdue and no other community seems to have as capable and flexible tools for it. Most of what I have is in these packages, but I've learnt quite a lot more since creating those:

https://github.com/mdsumner/graticule
https://github.com/mdsumner/tissot

More soon!