A more detailed look at spatial interpolation, integration through time, and I/O. For additional documentation e.g. see 1, 2, 3, 4. Here we illustrate a few things in more detail:
float_traj*data
)U,V
from gridded output to individual locationsu,v
from float_traj*data
OrdinaryDiffEq.jl
x(t),y(t)
from float_traj*data
using IndividualDisplacements, MITgcm
import IndividualDisplacements.OrdinaryDiffEq as OrdinaryDiffEq
import IndividualDisplacements.DataFrames as DataFrames
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/more/example123.jl"))
from MITgcm/pkg/flt
IndividualDisplacements.flt_example_download()
dirIn=IndividualDisplacements.flt_example_path
prec=Float32
df=read_flt(dirIn,prec);
using MeshArrays.jl
and e.g. a NamedTyple
𝑃,Γ=example2_setup();
plt=heatmap(Γ.mskW[1,1].*𝑃.u0,title="U at the start")
plt=heatmap(Γ.mskW[1,1].*𝑃.u1-𝑃.u0,title="U end - U start")
(select one trajectory)
tmp=df[df.ID .== 200, :]
tmp[1:4,:]
Super-impose trajectory over velocity field (first for u ...)
x=Γ.XG.f[1][:,1]
y=Γ.YC.f[1][1,:]
z=transpose(Γ.mskW[1].*𝑃.u0);
Super-impose trajectory over velocity field (... then for v)
x=Γ.XC.f[1][:,1]
y=Γ.YG.f[1][1,:]
z=transpose(Γ.mskW[1].*𝑃.v0);
dx=Γ.dx
uInit=[tmp[1,:lon];tmp[1,:lat]]./dx
nSteps=Int32(tmp[end,:time]/3600)-2
du=fill(0.0,2);
Visualize and compare with actual grid point values -- jumps on the tangential component are expected with linear scheme:
tmpu=fill(0.0,100)
tmpv=fill(0.0,100)
tmpx=fill(0.0,100)
for i=1:100
tmpx[i]=500.0 *i./dx
dxdt!(du,[tmpx[i];0.499./dx],𝑃,0.0)
tmpu[i]=du[1]
tmpv[i]=du[2]
end
And similarly in the other direction
tmpu=fill(0.0,100)
tmpv=fill(0.0,100)
tmpy=fill(0.0,100)
for i=1:100
tmpy[i]=500.0 *i./dx
dxdt!(du,[0.499./dx;tmpy[i]],𝑃,0.0)
tmpu[i]=du[1]
tmpv[i]=du[2]
end
Compare recomputed velocities with those from pkg/flt
nSteps=2998
tmpu=fill(0.0,nSteps); tmpv=fill(0.0,nSteps);
tmpx=fill(0.0,nSteps); tmpy=fill(0.0,nSteps);
refu=fill(0.0,nSteps); refv=fill(0.0,nSteps);
for i=1:nSteps
dxy_dt_replay(du,[tmp[i,:lon],tmp[i,:lat]],tmp,tmp[i,:time])
refu[i]=du[1]./dx
refv[i]=du[2]./dx
dxdt!(du,[tmp[i,:lon],tmp[i,:lat]]./dx,𝑃,Float64(tmp[i,:time]))
tmpu[i]=du[1]
tmpv[i]=du[2]
end
Solve through time using OrdinaryDiffEq.jl
with
dxdt!
is the function computing d(position)/dt
uInit
is the initial condition u @ tspan[1]
tspan
is the time interval𝑃
are parameters for dxdt!
Tsit5
is the time-stepping schemereltol
and abstol
are tolerance parameterstspan = (0.0,nSteps*3600.0)
#prob = OrdinaryDiffEq.ODEProblem(dxy_dt_replay,uInit,tspan,tmp)
prob = OrdinaryDiffEq.ODEProblem(dxdt!,uInit,tspan,𝑃)
sol = OrdinaryDiffEq.solve(prob,OrdinaryDiffEq.Tsit5(),reltol=1e-8,abstol=1e-8)
sol[1:4]
Compare recomputed trajectories with originals from MITgcm/pkg/flt
ref=transpose([tmp[1:nSteps,:lon] tmp[1:nSteps,:lat]])
maxLon=80*5.e3
maxLat=42*5.e3
#show(size(ref))
for i=1:nSteps-1
ref[1,i+1]-ref[1,i]>maxLon/2 ? ref[1,i+1:end]-=fill(maxLon,(nSteps-i)) : nothing
ref[1,i+1]-ref[1,i]<-maxLon/2 ? ref[1,i+1:end]+=fill(maxLon,(nSteps-i)) : nothing
ref[2,i+1]-ref[2,i]>maxLat/2 ? ref[2,i+1:end]-=fill(maxLat,(nSteps-i)) : nothing
ref[2,i+1]-ref[2,i]<-maxLat/2 ? ref[2,i+1:end]+=fill(maxLat,(nSteps-i)) : nothing
end
ref=ref./dx;
This notebook was generated using Literate.jl.