Improvements and example for linesplitting (#228)

* improve split method (code from haakon-e) * fix split methods * fix previous commit * add working example to docs for lon_0=-160 * update/fix test * Update docs/src/index.md Co-authored-by: Haakon Ludvig Langeland Ervik <[email protected]> * Update docs/src/index.md Co-authored-by: Haakon Ludvig Langeland Ervik <[email protected]> * Update docs/src/index.md Co-authored-by: Haakon Ludvig Langeland Ervik <[email protected]> * replace with AbstractVector as suggested * fix * Change the definition of `split` to make the output a bit cleaner * Make the core `split` code a bit clearer and more generic --------- Co-authored-by: Haakon Ludvig Langeland Ervik <[email protected]> Co-authored-by: Anshul Singhvi <[email protected]>
MakieOrg · Jul 10, 2024 · 4e39e7e · 4e39e7e
1 parent 8c4c8cc
commit 4e39e7e
Show file tree

Hide file tree

Showing 3 changed files with 86 additions and 36 deletions.
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -74,6 +74,38 @@ fig
 As you can see, the axis automatically transforms your input from the `source`
 CRS (default `"+proj=longlat +datum=WGS84"`) to the `dest` CRS.
 
+### Changing central longitude
+
+Be careful! Each data point is transformed individually.
+However, when using `surface` or `contour` plots this can lead to errors when the longitudinal dimension "wraps" around the planet.
+
+To fix this issue, the recommended approach is that you (1) change the central longitude of the map transformation (`dest`), and (2) `circshift` your data accordingly for `lons` and `field`.
+
+```@example MAIN
+function cshift(lons, field, lon_0)
+   shift = @. lons - lon_0 > 180
+   nn = sum(shift)
+   (circshift(lons - 360shift, nn), circshift(field, (nn, 0)))
+end
+```
+
+```@example MAIN
+lons = -180:180
+lats = -90:90
+field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]
+
+lon_0 = -160
+(lons_shift, field_shift) = cshift(lons, field, lon_0)
+```
+
+```@example MAIN
+fig = Figure()
+ax = GeoAxis(fig[1,1]; dest = "+proj=eqearth +lon_0=$(lon_0)")
+surface!(ax, lons_shift, lats, field_shift, colormap=:balance)
+lines!.(ax, GeoMakie.coastlines(ax), color=:black, overdraw = true)
+fig
+```
+
 You can also use quite a few other plot types and projections:
 ```@example quickstart
 fieldlons = -180:180; fieldlats = -90:90

diff --git a/src/linesplitting.jl b/src/linesplitting.jl
@@ -9,8 +9,6 @@
 #
 #This is needed to fix e.g. coast line displays when lon_0 is not 0 but cutting polygons at lon_0+-180.
 
-Base.split(tmp::Vector{<:LineString},ga::GeoAxis) = @lift(split(tmp,$(ga.dest)))
-
 """
     coastlines(ga::GeoAxis)
 Split coastline contours when ga.dest includes a "+lon_0" specification.
@@ -19,53 +17,73 @@ coastlines(ga::GeoAxis)=split(coastlines(),ga)
 
 module LineSplitting
 
-	using GeometryBasics: LineString
-	using Makie: Observable, @lift
+	import GeometryBasics
+	import GeoInterface as GI, GeometryOps as GO
+	import Makie: Observable, @lift, lift
+	import GeoMakie: GeoAxis
+
     # Since we're overriding Base.split, we must import it
 	import Base.split
-
-	function regroup(tmp::Vector)
-		coastlines_custom=LineString[]
-		println(typeof(coastlines_custom))
-		for ii in 1:length(tmp)
-			push!(coastlines_custom,tmp[ii][:]...)
-		end
-		coastlines_custom
+
+	###
+	function split(tmp::GeometryBasics.LineString, lon0::Real)
+		# lon1 is the "antimeridian" relative to the central longitude `lon0`
+		lon1 = lon0 < 0.0 ? lon0+180 : lon0-180 
+		# GeometryBasics handles line nodes as polygons.
+		linenodes = GeometryBasics.coordinates(tmp)  # get coordinates of line nodes
+		# Find nodes that are on either side of lon0
+		cond = GI.x.(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 (TODO: this needs to inject a point at the appropriate place)
+		split_coords = @. view((linenodes,), UnitRange(start_inds, end_inds))  # For each start-end pair, get those coords
+		# reconstruct lines from points
+		split_lines = GeometryBasics.MultiLineString(GeometryBasics.LineString.(split_coords))
 	end
-
-	function split(tmp::Vector{<:LineString}, lon0::Real)
-		[split(a,lon0) for a in tmp]
+
+	function split(tmp::AbstractVector{<:GeometryBasics.LineString}, lon0::Real)
+		[split(a, lon0) for a in tmp]
+	end
+
+	###
+	split(tmp::GeometryBasics.LineString,dest::Observable) = @lift(split(tmp, $(dest)))
+
+	function split(tmp::AbstractVector{<:GeometryBasics.LineString}, dest::Observable)
+		@lift([split(a, $(dest)) for a in tmp])
 	end
+
+	###
+	split(tmp::GeometryBasics.LineString, ax::GeoAxis) = split(tmp, ax.dest)
 
-	function split(tmp::LineString, lon0::Real)
-		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 !any(==(3), tmp2) # no value in tmp2 is equal to 3
-			[tmp]
-		else # some value in tmp2 is equal to 3
-			jj=[0;findall(tmp2.==3)...;np+1]
-			[LineString(tmp[jj[ii]+1:jj[ii+1]-1]) for ii in 1:length(jj)-1]
+	function split(tmp::AbstractVector{<:GeometryBasics.LineString}, ax::GeoAxis)
+		lift(ax.scene, ax.dest) do dest
+			[split(a, dest) for a in tmp]
 		end
 	end
-
-	split(tmp::Vector,dest::Observable) = @lift(split(tmp, $(dest)))
-	split(tmp::Observable,dest::Observable) = @lift(split($(tmp), $(dest)))
-	split(tmp::Observable,dest::String) = @lift(split($(tmp), (dest)))
-
-	function split(tmp::Vector{<:LineString},dest::String)
+
+	###
+	function split(tmp::GeometryBasics.LineString, dest::String)
 		if occursin("+lon_0",dest)
 			tmp1=split(dest)
 			tmp2=findall(occursin.(Ref("+lon_0"),tmp1))[1]
 			lon_0=parse(Float64,split(tmp1[tmp2],"=")[2])
-			regroup(split(tmp,lon_0))
+			split(tmp,lon_0)
 		else
 			tmp
 		end
 	end
 
+	function split(tmp::AbstractVector{<:GeometryBasics.LineString},dest::String)
+		[split(a,dest) for a in tmp]
+	end
+
+	###
+
+#	split(tmp::Vector,dest::Observable) = @lift(split(tmp, $(dest)))
+	split(tmp::Observable,dest::Observable) = @lift(split($(tmp), $(dest)))
+	split(tmp::Observable,dest::String) = @lift(split($(tmp), (dest)))
+
 end
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -32,6 +32,6 @@ Makie.set_theme!(Theme(
         @test split(GeoMakie.coastlines(),"+lon_0=-160") isa Vector
         ga = GeoAxis(Figure();dest = "+proj=wintri +lon_0=-160")
         @test GeoMakie.coastlines(ga) isa Observable
-        @test GeoMakie.coastlines(ga)[][1] isa GeometryBasics.LineString
+        @test GeoMakie.coastlines(ga)[] isa AbstractVector
     end
 end