Proj has this projection implemented for Africa (and only Africa), but it turns out that it has been used by others, including Lee for a map of the Pacific, and Snyder for a US map. The advantage of this map projection is that it can be used to make maps with much less distortion than many of the normal map projections.
Lee’s Pacific map is no longer in publication, and the equations for it don’t seem to be in proj or any software I know of. It looks like this:
So as a fun little project I’ve started implementing the map projection in the Julia interface to GMT. I’ll probably make a showcase example when it’s done.
Edit: It is possibly this PROJ projection: Lee Oblated Stereographic - though the picture on the web page looks a bit different.
I’d already started coding up the projection in Julia before I discovered this PROJ option. The maths is kind of cool; it uses complex numbers. The extremely rough Julia code is below. First (Lon, lat) are transformed to stereographic coordinates in the complex plane, defaulting to a centre of projection at 165W, 10S. Then these coordinates are “Oblated” using the function in Lee’s paper.
function st_lonlat_to_z(
lon::Float64,
lat::Float64,
lon0::Float64 = -165.,
lat0::Float64 = -10.,
)
x = (tand(lat) - tand(lat0)*cosd(lon-lon0)) /
(tand(lat0)*tand(lat) + secd(lat0)*secd(lat) + cosd(lon-lon0))
y = secd(lat0)*sind(lon-lon0) /
(tand(lat0)*tand(lat) + secd(lat0)*secd(lat) + cosd(lon-lon0))
return x + y*im
end
function oblate(
z::ComplexF64,
a::Float64 = 10.0,
K::Float64 = 1.442632,
C::Float64 = 0.070530,
D::Float64 = 0.049386,
)
return a*(K*z + (C-D*im)*z^3)
end
z = st_lonlat_to_z(-175., -21.)
println(z)
w = oblate(z)
println(w)
Yeah, it does look kind of cool. And as the paper says - the Mercator is one of the worst possible projections for the Pacific, especially at the northern and southern extremes.
Here is the link to a page describing an Oblated Stereographic projection Snyder made for the US: 26. Modified-Stereographic Conformal projections | Eu, Mircea Supposedly the scale varies by only +/- 2% over the whole map of the US, which is pretty impressive.
The beauty of Julia is that it handles complex numbers out of the box. I was originally going to code this up in Fortran for that same reason, but then thought it was an opportunity to learn the Julia interface to GMT.
Yep, your equations seem spherical.
The problem with making maps are not the projections themselves, that have all PROJ projs available via mapproject. The problem is the frames and the clipping function that each projection must have. I recently revisited this when implementing the the Spilhaus projection. It’s because of that then when trying pscoast, I get
julia> coast(R=(0,270,-90,80), proj="+proj=lee_os", shore=true, show=1)
pscoast [ERROR]: Bad case in gmt_map_clip_path (500)
I think I know why the equations from the paper give different answers than PROJ. The x and y coordinates seem to be for a map “rotated” through 90 degrees from what we’d normally expect. Here’s the picture from the paper:
Initially I thought they’d rotated the image simply to make better use of the page. Now I think this is actually how the projection was intended to look.
And here is the confirmation that both PROJ and the spherical equations in the paper agree (with Earth radius 6,370,997 m - the value PROJ uses):
And from the Julia code: -1.9708373127574434e6 - 2.9634533517437507e6im
The two answers agree (after the 90 degree rotation thing) to within about 0.1m. That’s good enough for me. And as an added bonus, if one forces the floats in the Julia to be 32-bit instead of 64-bit, the answers still agree to less than a metre.
Neither. It’s a non-specialized way of GMT saying: this projection is not supported
The issue is that all supported projections have associated functions that compute the frames and the clip path that delimits the valid plotting domain (think on a polar map where the clip path is a circle). This is what makes an addition of a new PROJ projection (that is not already supported by GMT too) so work/debug demanding. Each proj must have a certain signature and the corresponding functions exists, even when they are not needed (as is the case is we don’t want to plot any frame/annotations.
I’ve yet to peruse the examples of using PROJ linked to above. But my initial thought for doing the coastlines was to project the coastline for the whole world (or a sufficient amount to cover the Pacific) into Cartesian (x, y) coordinates and simply use -JX to handle the clipping and frames (of course annotations with be in projected metres, not latitude and longitude). I was thinking of simply adding the geographic labels as text labels at the locations they need to be.
One small thing I was wondering about is why PROJ is still about 0.1m different (on the sphere) than the answers i got from the equations in Julia. It’s obviously numerical precision and all that - probably in the trigonometric functions. But someone must have solved this problem because for really precise applications of mapping I think they sometimes want millimetre precision. Just curious, that’s all.
I was expecting that he example I linked above would do most of this job, it not all. () it was created to do the Best Rectangular example. But, I was wrong, Likely this projection can’t be applied to the whole Earth (seems latitude bounded)
Zooming in seems to show what we are after on the very center.
The 0.1 difference is strange indeed and I don’t believe taht it’s do to a rounding issue. Computations in C(++) and Julia are exact/equal to the last bit of floating point.
Regarding the 0.1m difference - I have a suspicion it might be to do with the implementation of the trigonometric functions. A bit of searching found mention that these functions are often less accurate than the double precision arithmetic can handle. And I also saw something about Julia porting its trigonometric functions from something called openlibm. I’ve asked at work to see if anyone knows - ultra precise computations must be a thing for applications like GPS (just a guess), so I suspect it’s a known and solved problem.
And regarding the domain of the projection - almost certainly it cannot be applied to the whole Earth (which is why I’ll probably need to limit the coastline to the Pacific area). Snyder mentions something about this when discussing his projection for the US:
A world map using the GS50 (50-State) projection is almost illegible, with meridians and parallels intertwined like wild vines.
The coastlinesproj function accepts to constrain the longitudes (the limits argument) but not the latitude (which made no sense for the purpose it was developed for). Maybe you can change that easily.
I’m currently doing a massive processing of the 1 arcsecond ALOS AW3D30 DEM to plot the topography. Honestly, SRTM would be just fine too, but I wanted to see if I could process large amounts of the Japanese data.
I simplified the OpenStreetMap coastline with gmt simplify. A back of the envelope calculation suggested one pixel in my final map will be about 1 km, so I simplified the coastline so it was no more than 500 m from the original full resolution coastline. This massively reduced the file size to about 1% of the full resolution.
So you projected the entire coastlines and then limited the domain of of the projected coordinates? What -R did you use?
(At this scale you might as well had used the GMT coastlines)
