Merge pull request #144 from gaelforget/v0p6p1a

Gulf of Mexico example, API fix, DriftersDataset

Project.toml

@@ -1,7 +1,7 @@
 name = "Drifters"
 uuid = "bd752fb7-3f37-44cb-a8fb-f461137b623f"
 authors = ["gaelforget <[email protected]>"]
-version = "0.6.0"
+version = "0.6.1"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"

ext/DriftersMakieExt.jl

@@ -1,10 +1,10 @@
 module DriftersMakieExt
 
 	using Makie, Drifters
-	import Drifters: InDiPlot, DataFrame, demo
+	import Drifters: DriftersDataset, DataFrame, demo, gcdist
 	import Makie: plot
 
-	function plot(x::InDiPlot)
+	function plot(x::DriftersDataset)
 		if !isempty(x.options)
 			o=x.options
 			if string(o.plot_type)=="simple_plot1"
@@ -15,6 +15,8 @@ module DriftersMakieExt
 				global_plot1(x.data[:I],x.data[:df])
 			elseif string(o.plot_type)=="plot_start_end"
 				plot_start_end(x.data[:I])
+			elseif string(o.plot_type)=="jcon_drifters"
+				plot_drifters_jcon(x.data.gdf;x.options...)
 			else
 				println("unknown option (b)")	
 			end
@@ -31,7 +33,7 @@ module DriftersMakieExt
 ```
 using Drifters, CairoMakie
 include("basics/random_flow_field.jl")
-x=InDiPlot( data=(I=I,ϕ=ϕ), options=(plot_type=:simple_plot1,) )
+x=DriftersDataset( data=(I=I,ϕ=ϕ), options=(plot_type=:simple_plot1,) )
 plot(x)
 ```
 """
@@ -65,7 +67,7 @@ end
 ```
 using Drifters, CairoMakie
 include("basics/solid_body_rotation.jl")
-x=InDiPlot( data=(I=I,), options=(plot_type=:simple_plot2,) )
+x=DriftersDataset( data=(I=I,), options=(plot_type=:simple_plot2,) )
 plot(x)
 ```
 """
@@ -115,7 +117,7 @@ Plot initial and final positions, superimposed on a globalmap of ocean depth log
 using Drifters, GLMakie
 include("worldwide/global_ocean_circulation.jl")
 
-x=InDiPlot( data=(I=𝐼,df=tmp_🔴,), options=(plot_type=:global_plot1,) )
+x=DriftersDataset( data=(I=𝐼,df=tmp_🔴,), options=(plot_type=:global_plot1,) )
 fig,tt=plot(x)
 fig
 
@@ -185,4 +187,60 @@ end
 return fig
 end
 
+## Drifters as plotted in JuliaCon Proceedings paper
+
+EarthRadius=6371e3 #in meters
+res=1/2 #resolution if test case
+lola(x,y)=(-100+x*res,17+y*res) #convert x/y to lon/lat
+
+"""
+    plot_drifters_jcon(gdf ; prefix="",pol=[],xlims=(-180.0,180.0),ylims=(-90.0,90.0),vmax=10.0)
+
+```
+include("LoopCurrent_replay.jl")
+LoopC=DriftersDataset( data=(gdf=gdf,), options=(plot_type="jcon_drifters",
+				prefix=prefix,xlims=(-98,-78),ylims=(18,31),pol=pol) )
+plot(LoopC)
+```
+"""
+function plot_drifters_jcon(gdf ; 	plot_type="jcon_drifters", prefix="",pol=[],
+									xlims=(-180.0,180.0),ylims=(-90.0,90.0),vmax=10.0)
+    fi00()=Figure(size=(900,600),fontsize=24, textcolor=:grey90)
+	fi0=with_theme(fi00,theme_dark()) 
+    ax0=Axis(fi0[1,1],xlabel="longitude",ylabel="latitude",
+        title=prefix*"surface drifter trajectories and speed (m/s)")
+    cr=(0,2.5); cm=:speed
+
+    m=ones(maximum([size(D,1) for D in gdf]))
+    for D in gdf
+		if in("longitude",names(D))
+			x=D.longitude
+			y=D.latitude
+		else
+			tmp=lola.(D.x,D.y)
+			x=[x[1] for x in tmp]
+			y=[x[2] for x in tmp]
+		end
+
+        if length(x) > 10
+            dt=(in("t",names(D)) ? diff(D.t)[1] : diff(D.time)[1].value/1000)
+            v=EarthRadius/dt*[gcdist(x[i],x[i+1],y[i],y[i+1]) for i in [1:length(x)-1 ; length(x)-1]]
+
+            m.=1.0
+            sum(v.>vmax)>0 ? m[findall(v.>vmax)].=NaN : nothing
+            n=1:length(v)
+            lines!(x.*m[n],y.*m[n],color=v.*m[n],colorrange=cr, colormap=cm)
+        end
+    end
+
+    xlims!(xlims...); ylims!(ylims...)
+    !isempty(pol) ? lines!(pol,color=:mediumpurple,linewidth=4) : nothing
+	Colorbar(fi0[1,2], colorrange=cr, colormap=cm, height = Relative(0.65))
+
+	fi0
+end
+
+##
+
+
 end

src/API.jl

@@ -64,34 +64,25 @@ struct uvwArrays{Ty} <: FlowFields
 end
 
 """
-    FlowFields(; u::Union{Array,Tuple}=[], v::Union{Array,Tuple}=[], w::Union{Array,Tuple}=[], 
-    period::Union{Array,Tuple}=[], gridtype::Symbol=:centered)
+    FlowFields(;    u::Union{Array,Tuple}=[], v::Union{Array,Tuple}=[], w::Union{Array,Tuple}=[], 
+                    period::Union{Array,Tuple}=[])
 
 Construct FlowFields data structure based on keywords.
 
 ```
 uC, vC, _ = SimpleFlowFields(16)
+to_C_grid!(uC,dims=1)
+to_C_grid!(uC,dims=2)
 F=FlowFields(u=uC,v=vC,period=(0,10.))
 ```
 """
 function FlowFields(; u::Union{Array,Tuple}=[], v::Union{Array,Tuple}=[], w::Union{Array,Tuple}=[], 
-    period::Union{Array,Tuple}=[], gridtype::Symbol=:centered)
+    period::Union{Array,Tuple}=[])
     (isa(u,Tuple)||length(u[:])==2) ? (u0=u[1]; u1=u[2]) : (u0=u; u1=u)
     (isa(v,Tuple)||length(v[:])==2) ? (v0=v[1]; v1=v[2]) : (v0=v; v1=v)
     (isa(w,Tuple)||length(w[:])==2) ? (w0=w[1]; w1=w[2]) : (w0=w; w1=w)
     if isempty(period)
         @warn "period needs to be defined"
-    else
-        if gridtype==:centered
-            to_C_grid!(u0,dims=1)
-            to_C_grid!(u1,dims=1)
-            to_C_grid!(v0,dims=2)
-            to_C_grid!(v1,dims=2)
-            if !isempty(w0)
-                to_C_grid!(w0,dims=3)
-                to_C_grid!(w1,dims=3)
-            end
-        end
     end
     if !isempty(u0) && !isempty(v0)
         if !isempty(w0)
@@ -104,6 +95,22 @@ function FlowFields(; u::Union{Array,Tuple}=[], v::Union{Array,Tuple}=[], w::Uni
     end
 end
 
+"""
+    to_C_grid!(x;dims=0)
+
+```
+if gridtype==:centered
+    to_C_grid!(u0,dims=1)
+    to_C_grid!(u1,dims=1)
+    to_C_grid!(v0,dims=2)
+    to_C_grid!(v1,dims=2)
+    if !isempty(w0)
+        to_C_grid!(w0,dims=3)
+        to_C_grid!(w1,dims=3)
+    end
+end
+```
+"""
 to_C_grid!(x;dims=0) = begin
     if (dims==1)&&(ndims(x)==2)
         x.=0.5*(circshift(x, (1,0))+x)
@@ -488,8 +495,15 @@ function gcdist(I::Individuals)
     :lon => first => :lo1,:lon => last => :lo2,
     :lat => first => :la1,:lat => last => :la2)
 
-    gcdist(lo1,lo2,la1,la2) = acos(sind(la1)*sind(la2)+cosd(la1)*cosd(la2)*cosd(lo1-lo2))
     tmp.gcd=[gcdist(tmp.lo1[i],tmp.lo2[i],tmp.la1[i],tmp.la2[i]) for i in 1:size(tmp,1)]
     return tmp
 end
 
+"""
+    gcdist(lo1,lo2,la1,la2)
+
+Great circle distance (gcd in radians) between two positions.
+"""
+gcdist(lo1,lo2,la1,la2) = acos(min(1.0,max(-1.0,sind(la1)*sind(la2)+cosd(la1)*cosd(la2)*cosd(lo1-lo2))))
+
+EarthRadius=6371e3 #in meters

src/Downloads.jl

@@ -41,6 +41,10 @@ function __init__datadeps()
         ["https://github.com/gaelforget/flt_example/archive/v1.0.tar.gz"],
         ["3306d37f5e70d0fa253a5b87b527f727891b84179d8411971bbc10f3164d5fde"],      
         post_fetch_method=unpackDV))
+    register(DataDep("Gulf_of_Mexico","Sample Drifter Dataset for example in JuliaConProceedings2023",
+        ["https://zenodo.org/records/14229203/files/Drifters_example.jld2"],
+        ["8c6ae8865821ac261316bfea3f31d1d362f210f806a85bfe5bb9775f3f91c2fa"],
+        ))
 end
 
 """
@@ -57,6 +61,8 @@ function getdata(nam::String)
             datadep"global_ocean_circulation_outputs"
         elseif nam=="flt_example"
             datadep"flt_example"
+        elseif nam=="Gulf_of_Mexico"
+            datadep"Gulf_of_Mexico"
         else
             println("unknown dataset")
         end

src/Drifters.jl

@@ -14,13 +14,14 @@ include("initial_positions.jl")
 include("Downloads.jl")
 include("example_ECCO.jl")
 include("example_OCCA.jl")
+include("example_GOM.jl")
 
 DiffEqBase.solve!(I::Individuals,args...)=∫!(I::Individuals,args...)
 DataFrames.groupby(I::Individuals,args...) = groupby(I.🔴,args...)
 DataFrames.DataFrame(I::Individuals) = I.🔴
 
 export Individuals, ∫!, solve!, DataFrame, groupby
-export FlowFields, convert_to_FlowFields
+export FlowFields, convert_to_FlowFields, to_C_grid!
 export uvwArrays, uvArrays, uvwMeshArrays, uvMeshArrays
 export dxdt!, dxy_dt_CyclicArray, dxy_dt_replay
 export postprocess_MeshArray, add_lonlat!, postprocess_xy, interp_to_xy
@@ -36,15 +37,18 @@ end
 #flt_example_path = datadeps.getdata("flt_example")
 #flt_example_download() = datadeps.getdata("flt_example")
 
-abstract type AbstractDriftersDiagnostic <: Any end
+abstract type AbstractDriftersDataset <: Any end
 
-Base.@kwdef struct InDiPlot <: AbstractDriftersDiagnostic
+Base.@kwdef struct DriftersDataset <: AbstractDriftersDataset
     path :: String = tempdir()
     name :: String = "unknown"
     options :: NamedTuple = NamedTuple()
     data :: NamedTuple = NamedTuple()
 end
 
-export InDiPlot
+export DriftersDataset
+
+#backward compatibility:
+InDiPlot=DriftersDataset; export InDiPlot
 
 end # module

src/compute.jl

@@ -27,7 +27,7 @@ I=Individuals(F,pos...)
 ∫!(I)
 
 using CairoMakie
-J=InDiPlot( data=(I=I,), options=(plot_type=:simple_plot2,) )
+J=DriftersDataset( data=(I=I,), options=(plot_type=:simple_plot2,) )
 f=plot(J)
 
 ref=[6.16, 7.23, 1.29, 1.0] # hide

src/example_GOM.jl

@@ -0,0 +1,13 @@
+
+import JLD2
+
+function Gulf_of_Mexico_setup()
+    file=joinpath(datadeps.getdata("Gulf_of_Mexico"),"Drifters_example.jld2")
+    jld=JLD2.jldopen(file)
+    pol=jld["polygons"]
+    u=jld["u"]; v=jld["v"]
+    x0=jld["x0"]; y0=jld["y0"]
+    dT=jld["dT"]; nt=jld["nt"]
+    T=(0.0,dT)
+    (x0,y0,u=u,v=v,T=T,dT=dT,nt=nt,pol=pol)
+end

test/runtests.jl

@@ -32,7 +32,7 @@ MeshArrays.GridLoad(MeshArrays.GridSpec("PeriodicChannel",MeshArrays.GRID_LL360)
     xlims=(-85.0,5.0)
     ylims=(20.0,67.0)
 
-    x=Drifters.InDiPlot( data=(I=I,df=tmp_🔴,), options=(plot_type=:global_plot1,) )
+    x=Drifters.DriftersDataset( data=(I=I,df=tmp_🔴,), options=(plot_type=:global_plot1,) )
     fig,tt=CairoMakie.plot(x)
     @test isa(fig,CairoMakie.Figure)
 end
@@ -47,7 +47,7 @@ end
 	T=(0.0,10*86400.0)
 	∫!(I,T)
 
-    fig=CairoMakie.plot( InDiPlot( data=(I=I,), options=(plot_type=:plot_start_end,) ) )
+    fig=CairoMakie.plot( DriftersDataset( data=(I=I,), options=(plot_type=:plot_start_end,) ) )
     @test isa(fig,CairoMakie.Figure)
 end
 
@@ -73,7 +73,7 @@ end
     I=Individuals(F,x,y)
     solve!(I,T)
 
-    fig=CairoMakie.plot( InDiPlot( data=(I=I,ϕ=ϕ), options=(plot_type=:simple_plot1,) ) )
+    fig=CairoMakie.plot( DriftersDataset( data=(I=I,ϕ=ϕ), options=(plot_type=:simple_plot1,) ) )
 
     @test isa(fig,CairoMakie.Figure)
 end
@@ -132,6 +132,9 @@ end
     I=(position=zeros(3,2),ID=1:2,record=deepcopy(df))
     I=Individuals(I)
     @test isa(I,Individuals)
+
+    GM=Drifters.Gulf_of_Mexico_setup()
+    @test isa(GM.T,Tuple)
 end
 
 @testset "doctests" begin