improvements and example for linesplitting

docs/src/index.md

@@ -73,21 +73,28 @@ fig
 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.
 
-E.g., if the data have the dimensions
+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
-lons = 0.5:359.5
-lats = -90:90
-field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats];
+cshift(lons, lats, field,lon_0)=begin
+   nn=sum((lons.-lon_0).>180)
+   ([(lons[end-nn+1:end].-360)...,lons[1:end-nn]...],
+   	lats,circshift(field,(nn,0)))
+end
 ```
-a `surface!` plot with the default arguments will lead to artifacts if the data along longitude 179 and 180 have significantly different values.
-To fix this, there are two approaches: (1) to change the central longitude of the map transformation, by changing the projection destination used like so:
 
-```julia
-ax = GeoAxis(fig[1,1]; dest = "+proj=eqearth +lon_0=180")
+```@example MAIN
+lon_0=-160
+(x,y,c)=cshift(lons, lats, field,lon_0)
 ```
 
-_or_ (2), circshift your data appropriately so that the central longitude you want coincides with the center of the longitude dimension of the data.
+```@example MAIN
+fig = Figure()
+ax = GeoAxis(fig[1,1]; dest = "+proj=eqearth +lon_0=$(lon_0)")
+surface!(ax, x, y, c, colormap=:balance)
+lines!.(ax, GeoMakie.coastlines(ax),color=:black,overdraw = true)
+fig
+```
 
 ### Countries loaded with GeoJSON
 ```@example MAIN

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,71 @@ coastlines(ga::GeoAxis)=split(coastlines(),ga)
 
 module LineSplitting
 
-	using GeometryBasics: LineString
-	using Makie: Observable, @lift
+	import GeometryBasics
+	import Makie: Observable, @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
+
+	getlon(p::GeometryBasics.Point) = p[1]
+
+	###
+	function split(tmp::GeometryBasics.LineString, lon0::Real)
+		lon0<0.0 ? lon1=lon0+180 : lon1=lon0-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
-	
-	function split(tmp::Vector{<:LineString}, lon0::Real)
+
+	function split(tmp::Vector{<:GeometryBasics.LineString}, lon0::Real)
 		[split(a,lon0) for a in tmp]
 	end
-
-	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]
-		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)))
+	###
+	split(tmp::GeometryBasics.LineString,dest::Observable) = @lift(split(tmp, $(dest)))
+
+	function split(tmp::Vector{<:GeometryBasics.LineString},dest::Observable)
+		@lift([split(a,$(dest)) for a in tmp])
+	end
 
-	function split(tmp::Vector{<:LineString},dest::String)
+	###
+	split(tmp::GeometryBasics.LineString,ax::GeoAxis) = split(tmp, ax.dest)
+
+	function split(tmp::Vector{<:GeometryBasics.LineString},ax::GeoAxis)
+		[split(a,ax.dest) for a in tmp]
+	end
+
+	###
+	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::Vector{<: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

test/runtests.jl

@@ -27,8 +27,8 @@ Makie.set_theme!(Theme(
     @testset "Line Splitting" begin
         @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 Vector
+        @test GeoMakie.coastlines(ga)[1] isa Observable
     end
 
     # @testset "Examples" begin