mapImage() should obey filledContour arg

[dankelley]

A related issue got closed 9 months ago. Reopening now, although I'm not sure how (or if) to handle this, without reinventing a wheel.

[richardsc]

Here's an idea. Presumably the internal "filled contour" functions (looks like .filled.contour(), and C_filledcontour) require a rectangular matrix, which is the problem with the map projections.

What if the projection transformed variables (i.e. the plot x and y) are simply gridded on a regular grid that is fine enough to smooth any funniness that may occur due to the lonlat-to-xy transformation? This might be a good use of the boxcarAverage function, as something like interpBarnes would likely be too slow (and is probably overkill in terms of an algorithm).

Not sure if I communicated this clearly or not, but it might be a start ...

[richardsc]

Here's a bit of an example, at least of the projected filled contours. Labels might be tricky, or might require a way to overplot the imagep on the same axes as had already been drawn by the mapPlot command (note that example below creates a new plot with the imagep call, which presumably has no associated projection information).

rm(list=ls())
library(oce)

data(coastlineWorld)
data(topoWorld)

lon <- expand.grid(topoWorld[['longitude']], topoWorld[['latitude']])[,1]
lat <- expand.grid(topoWorld[['longitude']], topoWorld[['latitude']])[,2]

mapPlot(coastlineWorld)

xy <- mapproject(lon, lat)

xrange <- par('usr')[1:2]
yrange <- par('usr')[3:4]
xg <- seq(xrange[1], xrange[2], 0.01)
yg <- seq(yrange[1], yrange[2], 0.01)

tg <- boxcarAverage2D(xy$x, xy$y, as.vector(topoWorld[['z']]), xg, yg)$average

imagep(xg, yg, tg, filledContour=TRUE, col=oceColorsJet, breaks=10)
mapPolygon(coastlineWorld)

Which makes:

[dankelley]

A worrisome aspect of akima is in the doc below (from http://www.inside-r.org/packages/cran/akima/docs/interp) which suggests the code is in fortran (meaning it is as fast as it can get, likely).

References Akima, H. (1978). A Method of Bivariate Interpolation and Smooth Surface Fitting for Irregularly Distributed Data Points. ACM Transactions on Mathematical Software 4, 148-164.

Akima, H. (1996). Algorithm 761: scattered-data surface fitting that has the accuracy of a cubic polynomial. ACM Transactions on Mathematical Software 22, 362--371.

Note interp is a wrapper for the two versions interp.old (it uses (almost) the same Fortran code from Akima 1978 as the S-Plus version) and interp.new (it is based on new Fortran code from Akima 1996). For linear interpolation the old version is choosen, but spline interpolation is done by the new version.

[richardsc]

Update: The akima interp() is not quite as slow as I thought. As with all of these things, it depends on how many points are being interpolated. I tried an example using some reanalysis data (ERA, 0.7 deg grid) over the Nordic Sea region (~3000 points to interpolate), and it works quite nicely. Notice that none of the "map axes" are there, because I'm just using imagep to plot it. Also, it's been stretched from the original map plot because imagep sets mar differently.

However, as a test case, using the topoWorld dataset (O(200000) points), even if only doing the interpolation of a small portion (i.e. the same Nordic Seas region), the time is unreasonably slow (I didn't bother timing, as it took long enough that I just kept killing the window after several minutes). One solution would be to trim out any values outside of the map range before interpolating, I suppose.

[richardsc]

Here's an example based on restricting the size of the input data for interpolation, with topoWorld. The whole thing takes about a second, which is reasonable.

[dankelley]

A couple of thoughts.

1. binMean and binApply could trim (x,y,f) data outside the "pretty" box, before doing any other processing. That will make it possible to work with larger length(x). Let n be the number of x data. The proposed addition will cost of order n, i.e. better than the present scheme with cost nN for an N by N grid. These costs are just for trimming, of course, which may be small compared to the overall cost. The overall cost of calculation will be reduced according to how many data are discarded, and this could be a huge gain if people use worldwide data to look in a box only a few hundred km wide.

2. I think the real solution to this problem is going to involve one of the two possibilities.

   2.1. Maybe the build-in code can be modified. I've not looked at it so I cannot guess whether this looks simple or hard. (Update: code is at ~/src/R-3.0.2/src/library/graphics/src/plot3d.c)

   2.2. An xy grid could be constructed by interpolation within polygons, then followed by normal filled contours. Before doing this, I'd want to think about how possible algorithms would scale. A simple way is to fit the 4 points within a polygon to a model f=a+b*x+c*y+d*xy. Another way involves breaking a rectangle up into 4 triangles from a constructed midpoint, then using f=a+b*x+c*y in those triangles. The former has some problems, and the latter is more standard, and is used commonly used to trace contours (e.g. in Gri and matlab, and I imagine in R).

[richardsc]

I don't know if this is useful for this issue or not, but here is a function I have been using to do filled contours on map projections. It uses the akima package for 2D interpolation. I'll add it to the oce-issues repo, too:

mapFilledContour <- function(lon, lat, f, breaks, n=500, contour=FALSE, clwd=1.25, ccol=1, col=oceColorsJet, ...)
{
    require(akima)
    fv <- as.vector(f)
    lo <- expand.grid(lon, lat)[,1]
    la <- expand.grid(lon, lat)[,2]
    good <- !is.na(fv)
    lo <- lo[good]
    la <- la[good]
    fv <- fv[good]
    fi <- interp(lo, la, fv, seq(min(lo), max(lo), length.out=n), seq(min(la), max(la), length.out=n))

    nbreaks <- length(breaks)
    col <- col(nbreaks-1)
    zlim <- range(breaks)
    levels <- seq(zlim[1], zlim[2], length.out=nbreaks)

    mapImage(fi, col=col, zlim=zlim)
    mapContour(fi$x, fi$y, fi$z, levels=breaks[-nbreaks][-1], col=col, lwd=2, ...)
    if (contour) mapContour(fi$x, fi$y, fi$z, levels=breaks[-nbreaks][-1], col=ccol, lwd=clwd, ...)
}

[richardsc]

Unfortunately it's ridiculously slow for largish data sets (such as topoWorld). It has been working well for contouring ERA data (0.7 degree) over the Nordic Sea region, but my first trial example using topoWorld has broken my computer ...

[richardsc]

While it's nice to close off issue to keep the list clean, I for one would argue that we should leave this open, or, at least start a new (cleaner) issue, because:

1. The issue has not been resolved, and as such currently oce spits a warning, and
2. In my view, this functionality is essential for oceanographic plotting capabilities in R to be on par with similar tools, including MatlabTM and Python.

As I currently "require" this functionality, I am using a work-around function based on some of the code from earlier in the thread, but it is far from ideal (i.e. very slow and not exactly correct).

Now, that's not to say that I think this should be something that we (or you) devote significant resources to immediately. Why not just let it stay, as a reminder, and when we have the chance perhaps we can poke at it some more. I had plans to learn better how the actually filled contour code works, so that it could perhaps be ported to projection space.

[dankelley]

I think this may work now. Below is test code from the help page on mapImage(), and the results. One thing I don't like too much: the akima license is raising red flags in the build process, and I hope that doesn't mean CRAN will complain.

Package has a FOSS license but eventually depends on the following
package which restricts use:
  akima

Test code

library(oce)
data(coastlineWorld)
data(topoWorld)

par(mfrow=c(2,1), mar=c(2, 2, 1, 1))
longitudelim <- c(-70,-50)
latitudelim <- c(40,50)
topo <- decimate(topoWorld, by=3) # coarse to illustrate filled contours
topo <- subset(topoWorld, latlim[1] < latitude & latitude < latlim[2])
topo <- subset(topo, (360+lonlim[1]) < longitude & longitude < (360+lonlim[2]))
mapPlot(coastlineWorld, type='l',
        longitudelim=longitudelim, latitudelim=latitudelim,
        proj="polyconic", orientation=c(90,-60,0), grid=TRUE)
breaks <- seq(-5000, 0, 500)
mapImage(topo, col=oceColorsJet, breaks=breaks)
mapLines(coastlineWorld)
box()
mapPlot(coastlineWorld, type='l',
        longitudelim=longitudelim, latitudelim=latitudelim,
        proj="polyconic", orientation=c(90,-60,0), grid=TRUE)
mapImage(topo, filledContour=TRUE, col=oceColorsJet, breaks=breaks)
box()
mapLines(coastlineWorld)

Results

[richardsc]

The license issue is a bit weird, seeing as you're loading a a package that is already on CRAN. The license (ACM) is here:

http://www.acm.org/publications/policies/softwarecrnotice

Looks pretty FOSS to me, but I'm not an expert.

[dankelley]

It seems weird to me too. The ACM license is permitted by CRAN, but the check engine issues a "NOTE" if you use a package that employs it. You're supposed to avoid NOTEs in packages.

However, it turns out to be a bit advantageous that I don't make Oce depend on it directly because this way, if you don't want filled map contours, you don't have to install the akima package. I would do the same for mapproj and proj4 except that I want them to be installed so the package-testing phase will test them.

[richardsc]

I've tested this a bit, and it seems to work really well. I especially like the filledContour=number option to control the grid refinement, as well as the nicely written warning in the docs.

[dankelley]

Thanks. I think I'll close this. It can always be reopened. The best would be a new issue, because the precise issue here (make something happen when the argument is supplied) has been addressed, so a new issue should be about something incorrect about how the argument is handled, or how the drawing is done, etc.

Thanks for the patience on this issue, @richardsc; 356 is a long step from the current issue pointer!