Plotting magnetic declination contours on ortho projection

[chrisredfern]

Hello, A couple of related questions:
(1) I'm really pleased to be able to plot magnetic declination contours on global maps! I would like to plot mag dec contours on orthographic plots centred on longitude 0 and 90 degrees.
This works nicely:
mapPlot(coastlineWorld, projection="+proj=ortho +lon_0=90", col="lightgray")
But at +lon_0=0 it's not quite so good on a couple of landmasses (next to Scandinavia, and Antarctica); is that issue soluble (or has been solved in latest oce code)?
(2) I can plot contours on a robinson projection with the code below, but the second line gives an error on the ortho projection:

mapContour(lon, lat, dec, col='blue', levels=seq(-180, -5, 5),

•        lty=3, drawlabels=FALSE)
  

mapContour(lon, lat, dec, col='blue', levels=seq(-180, -20, 20))
Error in if (dist > 1.4 * w && dx != 0) { :
missing value where TRUE/FALSE needed

Is this simple for me to correct?
Many thanks!
Chris

[dankelley]

Thanks for raising this. You are not alone! Map projections are difficult.

First, regarding whether more recent versions act differently -- the simplest answer is to try installing the github version and seeing for yourself. Or, provide a reproducible example so someone else can try it. (The first option is best because it will only take you 5 minutes to try the latest version.)

Second, regarding problems with projections: we have had these for years. It is a tough problem to crack because we want a general solution that will work across a range of projections. Some of the oldest issues in this project relate to projections. Since it is tricky, there is really very little chance that a repair can be done for your own particular case.

Third, regarding the error, please open an issue and provide a reproducible example, using the reprex package. That way, the developers can try to reproduce your bug. The core developers are both in very busy times lately, so we cannot promise a quick solution to the bug.

Finally, there is a way to get any projection to sort of work. That is, draw topography, instead of coastlines. Oce comes with a coarse topography, but you can always download a much finer dataset and use that. (There are hints about doing all this in the oce docs.)

Below is a reprex that shows a way to escape the projection-of-coastline problem: draw bathymetry instead! In my own field, this is more relevant anyway: national boundaries have no importance to oceanography, and the coastline is only slightly relevant: depth is what matters.

PS. I hope you can see how this is nice, because it shows code that certainly works, along with the output ... so please use a reprex if you make an issue, which I recommend. Also, in your issue, include the output from sessionInfo() so we can know more about your setup. If we need to change oce, you'll need to set up your system to build from source, which means having some compilers installed, etc.

library(oce)
#> Loading required package: gsw
data(coastlineWorld)
data(topoWorld) # download and read your own data for finer resolution
latHalifax <- 44.6
lonHalifax <- -63.6
par(mar=rep(1, 4))
proj <- paste0("+proj=ortho +lon_0=", lonHalifax, " +lat_0=", latHalifax)
mapPlot(coastlineWorld, proj=proj, type="n",
    col=2, pch=20, cex=0.5, axes=FALSE,
    grid=FALSE, longitudelim=lonHalifax+c(-180,180), latitudelim=c(45,135))
mapImage(topoWorld, colormap=colormap(name="gmt_globe"))
mapLines(coastlineWorld)
mapGrid(col="white", lwd=1.4, polarCircle=5)
mapGrid(col="magenta", lwd=1.4, polarCircle=5)
box()
mapPoints(lonHalifax, latHalifax, col=2, pch=20, cex=3)

^{Created on 2021-09-30 by the reprex package (v2.0.1)}

[dankelley]

Oh, sorry, you'd want to remove my mapLines() on the coastline because, depending on the exact view, they can be ruined. I only had them here because I was using this in a lecture about the Coriolis force, and I wanted to point out both Halifax and Toronto ... and I needed the country boundary to guide my students to the latter city.

[chrisredfern]

Dear Dan,
Very many thanks indeed for the quick reply! Lovely image and code!
I'll work through that, re-install the latest oce from Github and if I can't work out the magnetic declination mapContour issue will post a reprex and sessionInfo().
Chris

[chrisredfern]

Just want to thank you very much for your suggestions and code! It produced a beautiful plot (attached) to which I've added magnetic declination field lines (from your code elsewhere on the site). I can now add my data to this. This has really brought tears to my eyes!!
The mapLines have worked ok. The problem I had with the mapContour code I tracked down to the 'drawlabels' argument; when set to FALSE the code runs ok so presumably it is just projection maths issue with label orientation. I can add relevant ones manually at a later stage so it is not a problem.

mag-dec_plot

[dankelley]

Chris,

Oh, that's a lovely plot! Thanks for sharing it.

Dan.