filled contours in `mapImage()`

[richardsc]

Feature request, please! :)

[dankelley]

Can you make up a sample file? I think this won't be too hard but it would help to have a sample file.

[richardsc]

Not sure what you mean ... Perhaps the Levitus examples in the other issue (about mapImage)? Or even with a topo file. What I'm thinking of is something like:

library(oce)
data(coastlineMaritimes)
data(topoMaritimes)
mapPlot(coastlineMaritimes)
mapImage(topoMaritimes, filledContour=TRUE)

Much like the way imagep currently works if you want filled contours instead of an image plot.

[dankelley]

Yes, that's all I needed -- someone else to type in test code :-)

[dankelley]

Hm... this is difficult to do. In imagep() I am using .filled.contour() to do the work. But this requires gridded data i.e. a longitude (say) vector that has length equal to one dimension of the z matrix, etc. This is not the case for map projections.

Put simply, the (polygonal) elements to be coloured on the page are not in a rectangular grid.

One approach would be to interpolate from the polygons onto a grid of high resolution, and then to fill contours of that. So long as the grid was fine compared to the visual acuity of the viewer -- what Apple calls "retina resolution", that would be OK.

I looked into ways to speed up the drawing of polygons, but that code is very snarled ... I see comments from the developers about its being too complicated to maintain.

[dankelley]

mapImage() now has a new argument called filledContour. However, it doesn't do anything yet, except spit out a warning.

Also, I did some more thinking on the idea of setting up a grid, and a further problem with that is that it will require a function to go from page coordinates back to lat-lon coordinates, which is unavailable with the maptools package (I think). It is available in the proj4 package, but the latter has some problems with common use cases (e.g. it cannot handle world-wide maps very well in some cases).

I don't expect to make significant progress on this issue, sorry.

[richardsc]

Ah, well. I figured it would be hard. Harder than I'm capable of, which is why I directed it to you :)

[dankelley]

Yeah. Let's leave this open. I just had another look at the proj4 library, and it has an xy-to-lonlat ability, but the thing is that it spews errors on coastlineWorld so I don't see it as being a good choice.

Still there may be other tricks I can play, e.g. if a polygon is smaller than some-such size (relative to a window) I could just use points() on it. That might speed things up.

I think the real solution is for me to dig into the source code for mapproj and see how it hooks into the proj4 code (which is the source of all of this), and devise a patch that does the reverse function. Then I could suggest the patch to the people who maintain mapproj; if they don't want it, I could put this projection directly into oce. (It's been open source for like 15 years, if memory serves, and if it is GPL then R can take it.)

[dankelley]

I tried profiling, using

library(oce)
data(coastlineWorld)
mapPlot(coastlineWorld, type='l',
        latitudelim=c(40,50), longitudelim=c(-70,-50),
        proj="polyconic", orientation=c(90, -60,0), grid=TRUE)
data(topoMaritimes)
Rprof("mapimage.out")
mapImage(topoMaritimes, col=oceColorsGebco)
Rprof(NULL)

and then did

head(summaryRprof("mapimage.out"))

to find

               self.time self.pct total.time total.pct
"polygon"          17.14    23.69      21.20     29.31
"mapproject"       14.74    20.38      34.30     47.41
"mapImage"         11.50    15.90      72.34    100.00
".C"                6.64     9.18       6.64      9.18
"xy.coords"         4.06     5.61       4.06      5.61

which indicates that the projection is also taking a lot of the time, nearly as long (14.74 units) as the polygon plotting (17.14 units). So, even if the polygon part could be done in zero time, the speedup would only be by 27 percent.

Of course, this can be deceptive. Maybe if I (say) collect grid points along a longitude line and hand them in en mass to mapproject, things will speed up, and ditto for polygon.

[dankelley]

mapImage() in R/map.R has a double loop that computes xy<-mapproject(...). (This is at line 493 of the file as it exists at the moment of typing this.)

If there is overhead in calling the function, it can be reduced by expanding lon and lat out into matrices, then computing xy for the grid points.

Actually, that would be more sensible anyhow ... so I'll do that when I get a chance, and instrument the code to see if it speeds things up.

[dankelley]

The last update here was a long time ago. I wonder if this issue is still in play?

[richardsc]

Well, it is for me in the sense that I would really love to be able to plot filled contours on a map projection. This is standard in other popular analysis packages (e.g. python through basemap, and matlab though m_map), and I think would be a great addition to the mapping functions in oce.

I did come up with a crude way of doing this awhile ago, using the suggestion from your comment up above:

One approach would be to interpolate from the polygons onto a grid of high resolution, and then to fill contours of that. So long as the grid was fine compared to the visual acuity of the viewer -- what Apple calls "retina resolution", that would be OK.

As I recall, it more or less worked, except for two things:

1. It was very slow to draw polygons
2. There were some issues with properly aligning the colormap breaks and the contour lines

[richardsc]

This issue has been duplicated by #356, so I'm going to close this because discussion has continued there.