Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

plotting across the 180th meridian #126

Open
haakon-e opened this issue Sep 16, 2022 · 6 comments
Open

plotting across the 180th meridian #126

haakon-e opened this issue Sep 16, 2022 · 6 comments

Comments

@haakon-e
Copy link
Contributor

haakon-e commented Sep 16, 2022

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

image

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 (#7 (comment))

@haakon-e
Copy link
Contributor Author

haakon-e commented Sep 16, 2022

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

image

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

@gaelforget
Copy link
Contributor

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.

Screen Shot 2022-09-27 at 7 06 32 PM

@haakon-e
Copy link
Contributor Author

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)

image

@haakon-e
Copy link
Contributor Author

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)

image

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
Copy link

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:

image

(also relevant to MakieOrg/Makie.jl#742)

@asinghvi17
Copy link
Member

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants