How to map rectangular lat/lon region on non-rectangular projection with GeoMakie?

Specific Domains Geo
, ,
1

I’m trying to map a rectangular lat/lon polygon on a non-rectangular map projection using GeoMakie,

When I use four corner points to define the region, GeoMakie draws straight lines between them. The example below shows what I am getting. I would like the polygon to follow the curvature of the meridians instead.

image
image800×600 33.4 KB

using Makie, CairoMakie
using GeoMakie
using Makie, CairoMakie
using GeoMakie

fig = Figure()
ax = GeoAxis(fig[1, 1], title="poly = [(60, -30),  (180, -30),  (180, 30),  (60, 30)]")
poly = [(60, -30), (180, -30), (180, 30), (60, 30)]
poly!(ax, poly, color=:red)
fig

The default projection for GeoMakie is Equal Earth as I understand it, and the source default is Plate Caree, so I thought this would work. The corner points are obviously projected properly on the map.

What am I missing? Is there some way to accomplish what I want?

2

The thing is just that Makie currently doesn’t subsample data when it encounters nonlinear projections. So the vertices of the poly are correct but the lines between them are straight in screen space. If you made a polygon with more points between the corners, it would look better. I’m not sure yet what a good general solution to this could be without making it slow for larger datasets

3

From the standpoint of Makie this makes sense.

For GeoMakie, I was thinking the dest keyword would set to plot a geodetic polygon. It would be nice to have plotting in GeoMake obey something like a transform keyword in cartopy. But maybe that’s not optimal or simple.

Might look at the GeoMakie code to see how it draws the graticules on non-rectangular projections, which curve as expected. That seems essentially the same problem.

4

Does GeometryOps.jl help you?

5

If you plot a Box with LatLon coords from Meshes.jl with viz it will automatically subsample the box. If you send the box to a Proj transform it will also subsample.

You can call the viz recipe on any axis, including the GeoAxis above.

6

You can always subsample manually, which will give you the correct polygon:

using Makie, CairoMakie
using GeoMakie

import GeometryOps as GO, GeoInterface as GI

fig = Figure()
ax = GeoAxis(fig[1, 1], title="poly = [(60, -30),  (180, -30),  (180, 30),  (60, 30)]")
poly = GI.LineString([(60, -30), (180, -30), (180, 30), (60, 30)])
spoly = GO.segmentize(poly, max_distance = 1) # a point every degree of lon/lat
poly!(ax, spoly, color=:red)
fig

You can also use geodesic sampling:

poly!(ax, GO.segmentize(GO.GeodesicSegments(max_distance = 100_000), poly), color = (:orange, 0.5))
fig

I guess this sat in my drafts for a bit…

1 Like