mapImage(..., filledContour=TRUE) fails for antarctica

[richardsc]

Not sure what the details of this are, but I'm just dumping the results here for now to look into later:

library(oce)
data(coastlineWorld)
data(topoWorld)
mapPlot(coastlineWorld, proj='+proj=stere +lat_0=-90', longitudelim=c(-180, 180), latitudelim=c(-90, -50))
mapImage(topoWorld, colormap=colormap(name='gmt_relief'))
mapImage(topoWorld, colormap=colormap(name='gmt_relief'), filledContour=TRUE)

Error in seq.default(rx[1], rx[2], length.out = f * length(longitude)) : 
  'to' cannot be NA, NaN or infinite
> traceback()
4: stop("'to' cannot be NA, NaN or infinite")
3: seq.default(rx[1], rx[2], length.out = f * length(longitude))
2: seq(rx[1], rx[2], length.out = f * length(longitude))
1: mapImage(topoWorld, colormap = colormap(name = "gmt_relief"), 
       filledContour = 100)

[dankelley]

I'm going to take notes in the present comment, so I can see problem and notes without scrolling.

• made an issue testcode
• fixed reported problem but there is also a problem is in the akima::interp() call.
• made a new branch called issue721 to work on this

[dankelley]

I've got this working better on the test code, but it is too slow for interactive use because the call to akima::interp() takes of order a minute on this case. A bit of web searching suggests that others have found problems with interp() on large datasets. Probably I should explore methods other than akima::interp(), perhaps using oce::binMean2D() as a first pass, to thin out the data.

[dankelley]

@richardsc how much do you actually want filledContour? Is it just a matter of its being fast? I ask because I have determined (on this test case) that it will not be much faster than the default case (of drawing the non-rectangular polygons for the lon-lat boxes). The screenshot shows (code not yet in GH). The left panel shows that the filledContour method is only twice as fast as the default method, but the right panel shows that it's results are pretty poor in this case, with (a) open regions because of meridian divergence (something I can fix with very little extra calculation) but also (b) really quite coarse gridding. Doing this sort of interpolation with akima:interp() avoids the gappy problem but it is far too slow to be practical in large-area views.

Right now, my plan is to (a) write gap-filling code and (b) document very clearly that the filledContour=TRUE method is only for rough work. Actually, I think the argument is poorly named; if I were a user, I'd prefer to see something like method="accurate" (which will do polygons) and method="fast" (which does as filledContour=TRUE). I don't know how many people are using this code as is, and so I'm inclined to suggest changing the argument. I can always add a "..." to it, and then check to see if filledContour is in names(list(...)) and show a deprecation notice...

Below is the screenshot.

[dankelley]

I've added gap-filling code (branch issue721, commit 02929cc65b0a7461e0df212cf84c977cef18ef07), which lowers the aesthetic distraction although it is showing interpolants instead data, which I don't like. When I get a chance I'm going to see why there are big green splotches near AA; the grid cells are too small for that to result from bin averaging, so they are a mystery to me. Below is what I have now ... but note that I am really getting quite unkeen on the whole idea of filledContour (maybe @richardsc can convince me that it's worth keeping?)

[dankelley]

@richardsc -- I think this issue (as described in the title) can be closed. I opened a new one for the colour problem.

[richardsc]

Sorry for just opening this somewhat vague issue and then disappearing. To be clear, this wasn't a high priority plot, but rather something that I encountered while trying various projections and views. Sometimes I like to just try something that someone might do (like, someone who works on the Southern Ocean), even if it's not a part of my work, to see what will break.

I will close the issue, and move over to the discussion in #726.

@dankelley said: but note that I am really getting quite unkeen on the whole idea of filledContour (maybe @richardsc can convince me that it's worth keeping?)

For what it's worth, I think that having decent (and useable) filled contours for map projections is a "critical" feature for oce. Maybe not necessarily for making topo plots (though it can look nice), but also for plotting oceanic and atmospheric fields in projected space. Both the m_map package in matlab and basemap/cartopy in python can handle such plots (and in my experience with much less issue than oce), and it is a plotting feature that any experienced oceanographer would expect.

[dankelley]

Thanks for the notes on the need for filledContour=TRUE. I just now pushed a new version of the issue721 branch that I think settles the back-and-forth I've done over the past 3 days.