Line splitting for shifted reference longitude (#128)

The main user-facing function this introduces is `coastlines(ga)`, which returns an Observable of split coastlines.

* code to allow/fix coastlines with shifted reference longitude

* rename LineSplit as split and call within GeoAxis

* fix test

* mv LineSplitting to src/utils.jl

* improve test for LineSplitting

* Update test/runtests.jl

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

* Update src/utils.jl

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

* Update src/utils.jl

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

* Update src/utils.jl

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

* updates

* treat case of observable

* Move the LineSplitting module to a separate file

* More precise imports

* Cover all Observable cases for `split`

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

* Short-circuit condition for splitting

* Fix tests

* Move all code changes to linesplitting.jl

* reorder loading

---------

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

src/GeoMakie.jl

@@ -52,6 +52,7 @@ include("makie-axis.jl")
 
 # some basic recipes
 include("mesh_image.jl")
+include("linesplitting.jl")
 
 export GeoAxis, datalims, datalims!, automatic
 

src/linesplitting.jl

@@ -0,0 +1,71 @@
+
+############################################################
+#                                                          #
+#     Splitting Of LineString / coast lines                #
+#                                                          #
+############################################################
+
+#`LineSplitting` module defines a `split` method for vectors of `LineString` objects.
+#
+#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.
+"""
+coastlines(ga::GeoAxis)=split(coastlines(),ga)
+
+module LineSplitting
+
+	using GeometryBasics: LineString
+	using Makie: Observable, @lift
+    # 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
+	end
+
+	function split(tmp::Vector{<: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)))
+
+	function split(tmp::Vector{<: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))
+		else
+			tmp
+		end
+	end
+
+end

test/runtests.jl

@@ -24,6 +24,13 @@ Makie.set_theme!(Theme(
         @test GeoMakie.coastlines()[1] isa GeometryBasics.LineString
     end
 
+    @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
+    end
+
     # @testset "Examples" begin
     #     geomakie_path = dirname(dirname(pathof(GeoMakie)))
     #     examples = readdir(joinpath(geomakie_path, "examples"); join = true)