plotting across the 180th meridian

[haakon-e]

Hi,

I'm having some issues understanding how to plot a localized plot in the range, say, 160E to 160W, i.e. across the 180th meridian. Ideally, I'd like to be able to set a limit like lonlims=(160, -160) and let GeoMakie interpret that instruction correctly.

I've made a test script testing various approaches:

using GLMakie, GeoMakie
function add_plot(lonlims, subfig)
    try
        ga = GeoAxis(subfig, coastlines = true, dest = "+proj=eqc",
            lonlims=lonlims,
            title="lonlims = $lonlims",
        )
        image!(ga, -180..180, -90..90, rotr90(GeoMakie.earth()); interpolate = false)
    catch e
        ga = Axis(subfig, title="lonlims = $lonlims")
    end
end

fig = Figure()
add_plot((-180, -130), fig[1,1])  # partial plot that doesn't cross 180/-180 works fine.
add_plot((-181, -130), fig[2,1])  # plotting west of that longitude doesn't work with this approach
add_plot((179, -130), fig[3,1])   # setting bounds in the range (-180, 180) leads to an error: 
# Invalid x-limits as xlims[1] <= xlims[2] is not met for (179.0f0, -130.0f0).
add_plot((-180, 180), fig[1,2])  # plotting entire globe works fine
add_plot((0, 360), fig[2,2])      # plotting from (0,360) leads to weird vertical axis
add_plot((-130, 179), fig[3,2])  # setting bounds in (-180, 180), but flipping relative to fig[3,1] works, but is "opposite" of what I want.

fig

I see this post may be related to #122, but the code is a bit convoluted and hard to reproduce. Inspired by #118, I tried experimenting with +lon_0= -180, or 180, but then GeoMakie.earth() disappeared entirely, and/or coastlines produced weird horizontal lines. It also appears it's advised against anyways (https://github.com/MakieOrg/GeoMakie.jl/issues/7#issuecomment-581078328)

[haakon-e]

Looking closer at the source code, and the circshift tip in #7, I figured out a hacky way to do this, but it messes up the x-axis coordinates (unsurprisingly):

using GLMakie, GeoMakie
fig = Figure()
lonlims = (-40, 40)
ga = GeoAxis(fig[1,1], 
# coastlines = true, 
dest = "+proj=eqc",
    lonlims=lonlims,
    # latlims = (-40, 40),
    title="lonlims = $lonlims",
)
# image!(ga, -180..180, -90..90, rotr90(GeoMakie.earth()); interpolate = false)
image!(ga, -180..180, -90..90, circshift(rotr90(GeoMakie.earth()), (360, 0)); interpolate = false)

## Code adjusted from `GeoAxis` construction:
coast_line = Makie.convert_arguments(PointBased(), GeoMakie.coastlines())[1]
coastline_attributes = (;)
coastplot = lines!(ga, coast_line; color = :black, coastline_attributes...)
translate!(coastplot, -200, 0, 99)  # left half of coastlines (0, 180)

coastplot = lines!(ga, coast_line; color = :black, coastline_attributes...)
translate!(coastplot, 200, 0, 99)  # right half of coastlines (-180, 0)
fig

It's also unclear to me why ±200 is the appropriate translation limits (but that isn't really important...)

[gaelforget]

An alternative way to circshift could be to add duplicate views on the array on both sides + extend the longitude range to -540..540. This would also accommodate those who may want to specify longitudes in the 0-360 range.

[haakon-e]

It appears to me you didn't create that plot with GeoAxis. I think it's a good suggestion in general, and in fact I experimented a bit with making it work with coastlines too, which I'll get back to. For this, I am struggling getting the zoom to work as I would expect. Consider e.g.:

julia> using GeoMakie, GLMakie

julia> fig = Figure()

julia> ga = GeoAxis(fig[1, 1]; 
               lonlims=(-360, 0), coastlines = false,  # setting lonlims range >360 doesnt work at all, btw. Here I try to "shift" the view to be centered on the 180th meridian (experience so far by adjusting the `dest` projection has been futile...
               dest = "+proj=eqc",
           )

julia> e1 = rotr90(GeoMakie.earth());

julia> img = image!(-540..180, -90..90, [e1; e1]; interpolate = false)  # duplicate earth view "left"
Image{Tuple{Vector{Float32}, Vector{Float32}, Matrix{RGB{Float32}}}}

julia> xlims!(ga, -240, -140)

[haakon-e]

In principle, I think something similar could be done with coastlines. As a barebones code if someone wants to pick it up, here's simple code to "translate" the entire coastlines by some amount of longitude offset:

(the split function is from #128)

"Translate the coastline, which is just translating each element of the coastline array"
function translate_coastline(coastline::Vector{<:LineString}, offset::Real)
    translate_coastline_element.(coastline, offset)
end

"Translate a single element of the coastline, i.e. a LineString"
function translate_coastline_element(celem::LineString, offset::Real = -360.0)
    translated_coords = coordinates(celem) .+ Point(offset, 0)
    return LineString(translated_coords)
end

getlon(p::Point) = p[1]

"Split a coastline element if it crosses lon0 (thanks to gaelforget for initial work, see #128 for details"
function split_direct(tmp::LineString,lon0=-160.0)::Vector{<:LineString}
    linenodes = coordinates(tmp)  # get coordinates of line nodes
    # Find nodes that are on either side of lon0
    cond = getlon.(linenodes) .>= lon0
    # Find interval starts and ends
    end_cond = diff(cond)  # nonzero values denote ends of intervals
    end_inds = findall(!=(0), end_cond)
    start_inds = [firstindex(linenodes);  end_inds .+ 1]  # starts of intervals
    end_inds = [end_inds; lastindex(linenodes)]  # ends of intervals
    # do the splitting
    split_coords = view.(Ref(linenodes), UnitRange.(start_inds, end_inds))  # For each start-end pair, get those coords
    # reconstruct lines from points
    split_lines = LineString.(split_coords)
    return split_lines
end

"Split the entire coastline, which is just a vector of LineStrings. Return a vector of LineStrings"
function split_direct(ls_vec::Vector{<:LineString},lon0=-160.0)::Vector{<:LineString}
    vcat(split_direct.(ls_vec, lon0)...)
end

As an example, here I translate the coastlines by -180 degrees: (see geoaxis.jl:123 for reference code):

julia> using GeoMakie, GLMakie

julia> fig = Figure()

julia> ga = GeoAxis(fig[1, 1]; 
               coastlines = false,  # plot coastlines manually
               dest = "+proj=eqc",
           )

julia> coast_line = GeoMakie.coastlines();

julia> trans_cl = translate_coastline(coast_line, -180);

julia> left_cl = split_direct(trans_cl, -180);  # split along 180th meridian since we are still in (-180, 180) view.

julia> coastplot = lines!(ga, left_cl; color = :black)

With this, one can in principle do

julia> coast_line = GeoMakie.coastlines();

julia> trans_cl = translate_coastline(coast_line, -360);

julia> combined_coastline = [trans_cl; coast_line];

julia> coastplot = lines!(ga, combined_coastline; color = :black)

However, since GeoMakie / Proj wraps the coordinates, I haven't been able to actually do this in practice ...

[briochemc]

FWIW, maybe there is some useful code at uxarray? They seem to have the tooling to "split" those cells that cross the wrapping longitude. See, e.g., these plots from their docs:

(also relevant to https://github.com/MakieOrg/Makie.jl/issues/742)

[asinghvi17]

Huh, thanks for the link! Will check it out.