mapContour draws false straight contour lines

[cfranzke72]

Hello,

I want to plot a map in Mollweide projection with contour lines. The problem is that at the eastern and western ends where the lines should go to the other end (in a way behind the map), mapContour draws straight lines to the other end though the map. Hard to describe in words, see the attached PDF file. Is there a way to avoid this? I also provide a link to my netcdf file and a minimal R code example below.

Thanks, Christian

Plot: velpotential.2200_2300.ensmean.850.pdf

Netcdf file: https://www.dropbox.com/scl/fi/ufd7zw4fk7fxm8mxb25p3/velpot.2200_2300.850.ensmean.nc?rlkey=zfttnairhx1072tivian6mzsc&dl=0

My minimal R code example: library(ncdf4) library(oce) library(rcartocolor) library(rcolors)

data(coastlineWorld)

Tfile<-nc_open("velpot.2200_2300.850.ensmean.nc")

lon=Tfile$dim[[3]]$vals lat=Tfile$dim[[4]]$vals tim=Tfile$dim[[1]]$vals

NLON=Tfile$dim[[3]]$len NLAT=Tfile$dim[[4]]$len Ntim=Tfile$dim[[1]]$len

eof=ncvar_get(Tfile)

cols<-get_color("RdBu",n=201) cols2<-get_color("PuOr",n=201)

pdf("velpotential.2200_2300.ensmean.850.pdf") par(mfrow=c(1,1),mar=c(0,4,15,2)) cm<-colormap(zlim=c(-2e7,2e7),col=cols) drawPalette(pos=1,col=cols2,zlim=c(-2e7,2e7)) mapPlot(coastlineWorld,projection="+proj=moll",grid=FALSE,col="lightgray",drawBox=FALSE) mapImage(lon,lat,eof,col=cols2,zlim=c(-2e7,2e7)) mapContour(lon,lat,eof,lwd=1.5) mapGrid(col=gray(0.2)) mapPolygon(coastlineWorld[['longitude']], coastlineWorld[['latitude']]) dev.off()

[dankelley]

Hi, and thanks for the report.

I'm afraid that there is no way around this sort of behaviour. It has to do with the methods employed to do the transformation from geographic coordinates into x-y coordinates on the page.

There are lots of problems with projections that all come down to this same thing. In fact, the oldest issues in oce relate to map projections for coastline shapes. Although I've had some ideas for solving the problem, I have not worked them through to a solution. I call this the "edge of the earth" problem. My goal is to solve it for all the oce-supported projections, which means, basically, detecting where the edge of the earth is. My present idea is to do that by making a grid in x-y space and then trying to inverse-project into lon-lat space. For points past the edge of the projected earth, this results in an error in the underlying projection software, and so my idea is to make a logical matrix for the grid in x-y space. Then I can find a convex hull to map back to the lon-lat space that designates the edge of the earth. Once I have that, I can use some software to trim the elements (coastline shapes or contour lines, etc.) and display them. That's the idea, anyway. But not all projections do the same thing for "off-world" points. The underlying projection software keeps improving but some of the projections I initially wanted to support in oce were removed from oce because trying to inverse-project offworld points led to an infinite loop, or an R crash, etc. It's all a bit tricky.

Anyway, the above is just for some background, to explain why I cannot offer a fix for you. Projections are hard, basically. Personally, and I know this will offend some oceanographers (including oce co-authors!), I would just contour in lon-lat space, without using a projection. That has the disadvantage of shape and area distortion, but projections have similar problems, trading off one problem for another. And there is an advantage to simple lon-lat space: a user interested in the value at a given location can find that location easily, as for normal plots.

Having said that, if you do want projections, perhaps try using mapImage() with a colormap that has levels that correspond to where you would put the contour lines.

Um ... after writing the above, I tried your code and I also get those lines. Then I thought to try the following simpler example, and it works. So, the above comments may not be germane to your problem. I have multiple meetings this morning so I won't be able to look into your problem more closely. I am wondering whether the problem is that your longitudes are very close to the 0,360 limits, or something like that. Some testing may be required.

@richardsc has lots of experience with projections -- we've talked about the edge-of-world problem many times over the years -- so may have some ideas.

library(oce)
#> Loading required package: gsw
data(coastlineWorld)
mapPlot(coastlineWorld, projection = "+proj=moll")
data(levitus, package="ocedata")
mapContour(levitus$longitude, levitus$latitude, levitus$SST)

^{Created on 2024-06-07 with reprex v2.1.0}

[dankelley]

Here's an example code that might be a starting point if you want to try the stepped-colour scheme. Note that I also close the netcdf file, and am making a png instead of a pdf, so it can be seen in this comment. I removed some unused variables, and am using the cmocean package for colours. Here I am making the colour breaks manually, but of course you can do e.g. max(abs(eof)) to get a range, then build a sequence from that ... just be sure to have it centred, so that the colourscheme won't have the "middle" point that is not zero in your variable.

PS. the cmocean package has lots of nice colour schemes, and I encourage you to try it out, since a lot of thought has gone into making things work well (see multiple articles on colour).

library(ncdf4)
library(oce)
library(cmocean)
data(coastlineWorld)

Tfile <- nc_open("velpot.2200_2300.850.ensmean.nc")
lon <- Tfile$dim[[3]]$vals
lat <- Tfile$dim[[4]]$vals
eof <- ncvar_get(Tfile)
nc_close(Tfile)

if (!interactive()) {
    # pdf("velpotential.2200_2300.ensmean.850.pdf")
    png("velpotential.2200_2300.ensmean.850.png",
        unit = "in", height = 4, width = 7, res = 300,
        pointsize = 11
    )
}
breaks <- seq(-6, 6, 1) * 1e6 # by inspection
par(mar = c(3, 1, 1, 2))
drawPalette(breaks = breaks, col = cmocean("balance"), pos = 1)
par(mar = c(3, 1, 1, 1))
mapPlot(coastlineWorld, projection = "+proj=moll", grid = FALSE, col = "lightgray", drawBox = FALSE)
mapImage(lon, lat, eof, breaks = breaks, col = cmocean("balance"))
mapGrid(col = gray(0.2))
mapPolygon(coastlineWorld[["longitude"]], coastlineWorld[["latitude"]])
if (!interactive()) {
    dev.off()
}

[richardsc]

Hi @cfranzke72 ! Thanks for using oce 😄

@dankelley gave lots of good context about why this is hard, but also hinted at a solution -- if you reformat your matrices so that they go from -180 to 180 instead of 0 to 360 it works just fine:

library(oce)
library(ncdf4)
data(coastlineWorld)

nc <- nc_open('velpot.2200_2300.850.ensmean.nc')
lon <- ncvar_get(nc, 'lon')
lat <- ncvar_get(nc, 'lat')
vel <- ncvar_get(nc, 'velopot')

II <- lon > 180
velnew <- rbind(vel[II, ], vel[!II, ])
lonnew <- lon - 180

mapPlot(coastlineWorld, col='grey')
mapContour(lonnew, lat, velnew)

[dankelley]

That rbind() trick is elegant, @richardsc.

[dankelley]

@cfranzke72 I wonder if you might be able to look at #2218, in which I lay out a "solution". I had thought of trying to make mapContour() realize when this problem might crop up, but decided against that approach, as explained at #2218.

If you think the present issue has been addressed, please close it. Otherwise, please add more comments to help us in the development of the package. I ask this because the convention in this project is for reporters, not developers, to close issues.

Thanks for the report, and your patience!

Dan.

[cfranzke72]

Thanks for your help. I can now make nice plots.