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.

Sign up for GitHub

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

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

Merged
merged 22 commits into from
May 18, 2024
Merged

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

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
f78c996
code to allow/fix coastlines with shifted reference longitude
gaelforget Sep 23, 2022
4130c12
rename LineSplit as split and call within GeoAxis
gaelforget Sep 24, 2022
c34ae02
Merge branch 'master' into LineSplitting
gaelforget Sep 24, 2022
fb425c6
fix test
gaelforget Sep 24, 2022
ba06630
mv LineSplitting to src/utils.jl
gaelforget Sep 24, 2022
75ebde2
improve test for LineSplitting
gaelforget Sep 24, 2022
046475e
Merge branch 'master' into LineSplitting
gaelforget Dec 13, 2022
aef4fb1
Update test/runtests.jl
gaelforget Mar 1, 2023
4a7bc90
Update src/utils.jl
gaelforget Mar 1, 2023
664ac73
Update src/utils.jl
gaelforget Mar 1, 2023
eeeb829
Update src/utils.jl
gaelforget Mar 1, 2023
79ea295
Merge branch 'master' into LineSplitting_20240517
gaelforget May 17, 2024
da2c8ef
updates
gaelforget May 17, 2024
75161af
Merge branch 'LineSplitting_20240517' into LineSplitting
gaelforget May 17, 2024
835686b
treat case of observable
gaelforget May 17, 2024
afe64ee
Move the LineSplitting module to a separate file
asinghvi17 May 18, 2024
aa8dbe9
More precise imports
asinghvi17 May 18, 2024
8868271
Cover all Observable cases for `split`
asinghvi17 May 18, 2024
9898014
Short-circuit condition for splitting
asinghvi17 May 18, 2024
c556885
Fix tests
asinghvi17 May 18, 2024
e2f7726
Move all code changes to linesplitting.jl
asinghvi17 May 18, 2024
4917e2e
reorder loading
asinghvi17 May 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/geoaxis.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,7 @@ function GeoAxis(args...;
Makie.Observables.connect!(ax.scene.transformation.transform_func, _transformation)

# Plot coastlines
coast_line = GeoMakie.coastlines()
coast_line = split(GeoMakie.coastlines(),dest)
coastplot = lines!(ax, coast_line; color = :black, coastline_attributes...)
translate!(coastplot, 0, 0, 99) # ensure they are on top of other plotted elements
xprot = ax.xaxis.protrusion[]
Expand Down
60 changes: 60 additions & 0 deletions src/utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -415,3 +415,63 @@ macro hijack_observable(name)
end)

end

############################################################
# #
# 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.

module LineSplitting
gaelforget marked this conversation as resolved.
Show resolved Hide resolved

import GeoMakie.LineString
import GeoMakie.Observable
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,lon0=-160.0)
gaelforget marked this conversation as resolved.
Show resolved Hide resolved
[split(a,lon0) for a in tmp]
end

function split(tmp::LineString,lon0=-160.0)
gaelforget marked this conversation as resolved.
Show resolved Hide resolved
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
Copy link
Contributor

@haakon-e haakon-e Nov 17, 2022 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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


split(tmp::Vector,dest::Observable) = tmp
gaelforget marked this conversation as resolved.
Show resolved Hide resolved

function split(tmp::Vector,dest::String)
gaelforget marked this conversation as resolved.
Show resolved Hide resolved
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
5 changes: 5 additions & 0 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,11 @@ Makie.set_theme!(Theme(
@test GeoMakie.coastlines()[1] isa GeometryBasics.LineString
end

@testset "LineSplitting" begin
gaelforget marked this conversation as resolved.
Show resolved Hide resolved
@test split(GeoMakie.coastlines(),"+lon_0=-160") isa Vector
@test split(GeoMakie.coastlines(),"+lon_0=-160")[1] isa GeometryBasics.LineString
end
gaelforget marked this conversation as resolved.
Show resolved Hide resolved

@testset "Examples" begin
geomakie_path = dirname(dirname(pathof(GeoMakie)))
examples = readdir(joinpath(geomakie_path, "examples"); join = true)
Expand Down