code to allow/fix coastlines with shifted reference longitude #128

gaelforget · 2022-09-23T00:44:34Z

This adresses

The new split method let's us go from

to

The new feature can be tested with:

function demo_GeoMakie(lon0=-160.0,fixIssue=true)
	f = Figure()
	ax = GeoMakie.GeoAxis(f[1,1]; dest = "+proj=longlat +datum=WGS84 +lon_0=$(lon0)")

	all_lines=split(GeoMakie.coastlines(),lon0)
	Antarctica=split(GeoMakie.coastlines()[99],lon0)

	if !fixIssue
		lines!(ax, GeoMakie.coastlines())
	else
		[lines!(ax, l,color=:black) for l in all_lines]
		lines!(ax, Antarctica[1],color=:green)
		lines!(ax, Antarctica[2],color=:red)
	end
	
	f
end

gaelforget · 2022-09-24T12:29:32Z

This should work with other projections such as

fig = Figure()
ax = GeoAxis(fig[1,1]; dest = "+proj=wintri  +lon_0=-160", coastlines=true)
fig

Note : issues remain with the axes as highlighted in this plot but that's separate from the coastline part I think.

gaelforget · 2022-11-04T15:35:28Z

Could I get a review on this PR?

@asinghvi17 maybe based on #118 (comment)

haakon-e · 2022-11-17T01:07:30Z

It seems GUI zooming doesn't quite work as expected, though I'm not entirely sure if that is actually related to this PR.

edit: after reviewing the PR more closely, it is clear that this is a separate issue.

If I set lon0=0, things work as expected:

however, setting lon0=160:

Note that I'm also on GLMakie#master, Makie#master, etc... but I can test this again on latest release.

haakon-e

Nice PR! I'm no expert on GeoMakie / Makie at all, but thought I'd leave some suggestions since the PR addresses some issues I had some time back.

haakon-e · 2022-11-17T03:46:34Z

src/utils.jl

+
+	function split(tmp::LineString,lon0=-160.0)
+		lon0<0.0 ? lon1=lon0+180 : lon1=lon0-180 
+		np=length(tmp)
+		tmp2=fill(0,np)
+		for p in 1:np
+			tmp1=tmp[p]
+			tmp2[p]=maximum( [(tmp1[1][1]<=lon1)+2*(tmp1[2][1]>=lon1) , (tmp1[2][1]<=lon1)+2*(tmp1[1][1]>=lon1)] )
+		end
+		if sum(tmp2.==3)==0
+			[tmp]
+		else
+			jj=[0;findall(tmp2.==3)...;np+1]
+			[LineString(tmp[jj[ii]+1:jj[ii+1]-1]) for ii in 1:length(jj)-1]
+		end
+	end


I think this can be improved syntactically (if I understand your code correctly), by:

Suggested change

function split(tmp::LineString,lon0=-160.0)

lon0<0.0 ? lon1=lon0+180 : lon1=lon0-180

np=length(tmp)

tmp2=fill(0,np)

for p in 1:np

tmp1=tmp[p]

tmp2[p]=maximum( [(tmp1[1][1]<=lon1)+2*(tmp1[2][1]>=lon1) , (tmp1[2][1]<=lon1)+2*(tmp1[1][1]>=lon1)] )

end

if sum(tmp2.==3)==0

[tmp]

else

jj=[0;findall(tmp2.==3)...;np+1]

[LineString(tmp[jj[ii]+1:jj[ii+1]-1]) for ii in 1:length(jj)-1]

end

end

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

function split(tmp::LineString,lon0=-160.0)

lon1 = lon0 + (lon0 < 0 ? 180 : -180)

linenodes = GeometryBasics.coordinates(tmp) # get coordinates of line nodes

# Find nodes that are on either side of lon0

cond = getlon.(linenodes) .>= lon1

# 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 = GeometryBasics.LineString.(split_coords)

end

My expectation is that concatenating all the split nodes should be the same as the original nodes, which this suggestion satisfies, i.e.

vcat(split_coords...) == linenodes

In your code however, it seems that if the input is

julia> tmp = GeoMakie.coastlines()[1] # longitudes of ≈ - 161 ± 2 julia> lon0 = 20 # i.e. I want to split the coastline along -160 julia> split_lines = split(tmp, lon0) julia> vcat(coordinates.(split_lines)) 16-element Vector{GeometryBasics.Point{2, Float32}}: [-163.71289, -78.595665] [-163.1058, -78.223335] [-163.1058, -78.223335] [-161.24512, -78.38018] [-161.24512, -78.38018] [-160.2462, -78.69365] [-159.4824, -79.04634] [-159.20819, -79.49701] [-161.1276, -79.63421] [-162.43985, -79.28146] [-162.43985, -79.28146] [-163.0274, -78.92877] [-163.0274, -78.92877] [-163.0666, -78.869965] [-163.0666, -78.869965] [-163.71289, -78.595665]

then there are many duplicate coordinates in the output.

Thank you very much for the suggested improvements. I tested your revised code back in Nov and it worked well from what I remember. Wanted to take another look to confirm I understood your code but then I got side tracked -- apologies for that.

Issue now though is that my basic test case no longer works (ax = GeoMakie.GeoAxis(f[1,1]; dest = "+proj=longlat +datum=WGS84 +lon_0=-160",lonlims=GeoMakie.automatic)). I assume it's because of code merged ahead of this one but I am not sure.

Have now merged master into branch and updated example code, so I could try suggested change again. Might wait to hear from @asinghvi17 though

src/utils.jl

test/runtests.jl

asinghvi17 · 2023-01-12T11:57:48Z

Will review this in a couple days.

asinghvi17

Oops - looks like I didn't post my review when I last looked at this. I'll look at it again in the next couple days.

Couple questions:

Why is this in a Module?
does this only work on LineStrings for now, or is there a more generic way?

Simon and I just spoke last night and we're trying to get this kind of splitting directly in the transform, so that we can override the behaviour for individual transforms as well (e.g., some of the weird transforms which aren't "convex").

src/utils.jl

test/runtests.jl

src/utils.jl

Co-authored-by: Anshul Singhvi <[email protected]>

Co-authored-by: Haakon Ludvig Langeland Ervik <[email protected]>

src/utils.jl

gaelforget · 2024-05-17T14:27:30Z

Oops - looks like I didn't post my review when I last looked at this. I'll look at it again in the next couple days.

Couple questions:

Why is this in a Module?

does this only work on LineStrings for now, or is there a more generic way?

Simon and I just spoke last night and we're trying to get this kind of splitting directly in the transform, so that we can override the behaviour for individual transforms as well (e.g., some of the weird transforms which aren't "convex").

Could we merge this (freshly updated) PR and punt further improvements to later revisions (when new issue arises)?

The current implementation is not that intrusive (see module), PR provides test & docstrings, and having this feature inside GeoMakie would make the package more potentially useful for ocean viz.

If a general solution were eventually found and implemented (for now this PR only treats the lon_0 & LineStrings case) then that's great, but this might take a few more years, who knows ...

Currently I am just carrying what I need in an extension to MeshArrays ( https://github.com/JuliaClimate/MeshArrays.jl/blob/master/ext/MeshArraysMakieExt.jl ), because I couldn't get this merged here in timely manner. That's less than ideal, so I am giving GeoMakie one more try now.

asinghvi17 · 2024-05-18T05:21:12Z

Sure, we can merge for now but I will refactor a bit to make it easier to merge #207. That PR also adds a dependency on GeometryOps.jl that should make it a lot easier to extend this for arbitrary geometry.

This should definitely be changed to an internal function at some point, though

gaelforget · 2024-05-18T11:22:44Z

Sure, we can merge for now but I will refactor a bit to make it easier to merge #207. That PR also adds a dependency on GeometryOps.jl that should make it a lot easier to extend this for arbitrary geometry.

Great! Thank you very much @asinghvi17 & @haakon-e

asinghvi17 · 2024-12-14T16:58:39Z

It looks like this doesn't do exactly what I had expected at first:

julia> ls = LineString(Point2f[(-179, 0), (179, 0)])
LineString{2, Float32}(Point{2, Float32}[[-179.0, 0.0], [179.0, 0.0]])

julia> split(ls, 0)
MultiLineString{2, Float32}(LineString{2, Float32}[LineString{2, Float32}(Point{2, Float32}[[-179.0, 0.0], [179.0, 0.0]])])

julia> split(ls, 180)
MultiLineString{2, Float32}(LineString{2, Float32}[LineString{2, Float32}(Point{2, Float32}[[-179.0, 0.0]]), LineString{2, Float32}(Point{2, Float32}[[179.0, 0.0]])])

so will refactor it using a different approach. (Did we miss a +180 to the longitudes somewhere @haakon-e? That seems to be required to make this work...)

haakon-e · 2024-12-14T20:54:02Z

@asinghvi17 It's a long time since I've thought about this PR, could you explain in a bit more detail what kind of output you would expect from each case above?

code to allow/fix coastlines with shifted reference longitude

f78c996

gaelforget mentioned this pull request Sep 23, 2022

lonlims seems to misbehave #118

Open

rename LineSplit as split and call within GeoAxis

4130c12

Merge branch 'master' into LineSplitting

c34ae02

gaelforget marked this pull request as draft September 24, 2022 12:43

gaelforget added 3 commits September 24, 2022 09:31

fix test

fb425c6

mv LineSplitting to src/utils.jl

ba06630

improve test for LineSplitting

75ebde2

gaelforget marked this pull request as ready for review September 24, 2022 13:54

haakon-e reviewed Nov 17, 2022

View reviewed changes

This was referenced Nov 17, 2022

GLMakie resizing is weird with lon0 != 0 #136

Open

plotting across the 180th meridian #126

Open

Merge branch 'master' into LineSplitting

046475e

asinghvi17 reviewed Feb 1, 2023

View reviewed changes

src/utils.jl Outdated Show resolved Hide resolved

test/runtests.jl Outdated Show resolved Hide resolved

src/utils.jl Outdated Show resolved Hide resolved

SimonDanisch mentioned this pull request Feb 1, 2023

turn GeoAxis into a block and fix projection #148

Closed

5 tasks

gaelforget and others added 7 commits March 1, 2023 15:36

Update test/runtests.jl

aef4fb1

Co-authored-by: Anshul Singhvi <[email protected]>

Update src/utils.jl

4a7bc90

Co-authored-by: Haakon Ludvig Langeland Ervik <[email protected]>

Update src/utils.jl

664ac73

Co-authored-by: Haakon Ludvig Langeland Ervik <[email protected]>

Update src/utils.jl

eeeb829

Co-authored-by: Haakon Ludvig Langeland Ervik <[email protected]>

Merge branch 'master' into LineSplitting_20240517

79ea295

updates

da2c8ef

Merge branch 'LineSplitting_20240517' into LineSplitting

75161af

asinghvi17 reviewed May 17, 2024

View reviewed changes

src/utils.jl Outdated Show resolved Hide resolved

treat case of observable

835686b

asinghvi17 added 7 commits May 18, 2024 10:56

Move the LineSplitting module to a separate file

afe64ee

More precise imports

aa8dbe9

Cover all Observable cases for split

8868271

This should definitely be changed to an internal function at some point, though

Short-circuit condition for splitting

9898014

Fix tests

c556885

Move all code changes to linesplitting.jl

e2f7726

reorder loading

4917e2e

asinghvi17 merged commit 0fd8350 into MakieOrg:master May 18, 2024
2 checks passed

gaelforget mentioned this pull request May 20, 2024

improvements and example for linesplitting #228

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

code to allow/fix coastlines with shifted reference longitude #128

code to allow/fix coastlines with shifted reference longitude #128

gaelforget commented Sep 23, 2022 •

edited

Loading

gaelforget commented Sep 24, 2022

gaelforget commented Nov 4, 2022

haakon-e commented Nov 17, 2022 •

edited

Loading

haakon-e left a comment

haakon-e Nov 17, 2022 •

edited

Loading

gaelforget Jan 13, 2023

gaelforget May 17, 2024

asinghvi17 commented Jan 12, 2023

asinghvi17 left a comment

gaelforget commented May 17, 2024 •

edited

Loading

asinghvi17 commented May 18, 2024

gaelforget commented May 18, 2024

asinghvi17 commented Dec 14, 2024

haakon-e commented Dec 14, 2024

-	function split(tmp::LineString,lon0=-160.0)
-		lon0<0.0 ? lon1=lon0+180 : lon1=lon0-180
-		np=length(tmp)
-		tmp2=fill(0,np)
-		for p in 1:np
-			tmp1=tmp[p]
-			tmp2[p]=maximum( [(tmp1[1][1]<=lon1)+2*(tmp1[2][1]>=lon1) , (tmp1[2][1]<=lon1)+2*(tmp1[1][1]>=lon1)] )
-		end
-		if sum(tmp2.==3)==0
-			[tmp]
-		else
-			jj=[0;findall(tmp2.==3)...;np+1]
-			[LineString(tmp[jj[ii]+1:jj[ii+1]-1]) for ii in 1:length(jj)-1]
-		end
-	end
+        getlon(p::GeometryBasics.Point) = p[1]
+	function split(tmp::LineString,lon0=-160.0)
+		lon1 = lon0 + (lon0 < 0 ? 180 : -180)
+		linenodes = GeometryBasics.coordinates(tmp)  # get coordinates of line nodes
+		# Find nodes that are on either side of lon0
+		cond = getlon.(linenodes) .>= lon1
+		# 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 = GeometryBasics.LineString.(split_coords)
+	end

code to allow/fix coastlines with shifted reference longitude #128

code to allow/fix coastlines with shifted reference longitude #128

Conversation

gaelforget commented Sep 23, 2022 • edited Loading

gaelforget commented Sep 24, 2022

gaelforget commented Nov 4, 2022

haakon-e commented Nov 17, 2022 • edited Loading

haakon-e left a comment

Choose a reason for hiding this comment

haakon-e Nov 17, 2022 • edited Loading

Choose a reason for hiding this comment

gaelforget Jan 13, 2023

Choose a reason for hiding this comment

gaelforget May 17, 2024

Choose a reason for hiding this comment

asinghvi17 commented Jan 12, 2023

asinghvi17 left a comment

Choose a reason for hiding this comment

gaelforget commented May 17, 2024 • edited Loading

asinghvi17 commented May 18, 2024

gaelforget commented May 18, 2024

asinghvi17 commented Dec 14, 2024

haakon-e commented Dec 14, 2024

gaelforget commented Sep 23, 2022 •

edited

Loading

haakon-e commented Nov 17, 2022 •

edited

Loading

haakon-e Nov 17, 2022 •

edited

Loading

gaelforget commented May 17, 2024 •

edited

Loading