Merge pull request #75 from JuliaClimate/v0p3p1

V0p3p1

Project.toml

@@ -1,7 +1,7 @@
 name = "IndividualDisplacements"
 uuid = "b92f0c32-5b7e-11e9-1d7b-238b2da8b0e6"
 authors = ["gaelforget <[email protected]>"]
-version = "0.3.0"
+version = "0.3.1"
 
 [deps]
 CFTime = "179af706-886a-5703-950a-314cd64e0468"
@@ -22,9 +22,9 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 CFTime = "0.1"
 CSV = "0.6, 0.7, 0.8"
 CyclicArrays = "0.2, 0.3, 0.4, 0.5"
-DataFrames = "0.21, 0.22"
-MITgcmTools = "0.1.6"
-MeshArrays = "0.2.14"
+DataFrames = "0.21, 0.22, 1.0, 1.1"
+MITgcmTools = "0.1.20"
+MeshArrays = "0.2.19"
 NetCDF = "0.10, 0.11"
 OceanStateEstimation = "0.1.5"
 OrdinaryDiffEq = "5"

examples/basics/detailed_look.jl

@@ -46,8 +46,8 @@ df=read_flt(dirIn,prec);
 # ## 4. Visualize Velocity Fields
 #
 # ```
-# plt=heatmap(Γ["mskW"][1,1].*𝑃.u0,title="U at the start")
-# plt=heatmap(Γ["mskW"][1,1].*𝑃.u1-𝑃.u0,title="U end - U start")
+# plt=heatmap(Γ.mskW[1,1].*𝑃.u0,title="U at the start")
+# plt=heatmap(Γ.mskW[1,1].*𝑃.u1-𝑃.u0,title="U end - U start")
 # ```
 
 # ## 5. Visualize Trajectories
@@ -59,25 +59,25 @@ 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);
+x=Γ.XG.f[1][:,1]
+y=Γ.YC.f[1][1,:]
+z=transpose(Γ.mskW[1].*𝑃.u0);
 
 # plt=contourf(x,y,z,c=:delta)
 # plot!(tmp[:,:lon],tmp[:,:lat],c=:red,w=4,leg=false)
 
 # Super-impose trajectory over velocity field (... then for v)
 
-x=Γ["XC"].f[1][:,1]
-y=Γ["YG"].f[1][1,:]
-z=transpose(Γ["mskW"][1].*𝑃.v0);
+x=Γ.XC.f[1][:,1]
+y=Γ.YG.f[1][1,:]
+z=transpose(Γ.mskW[1].*𝑃.v0);
 
 # plt=contourf(x,y,z,c=:delta)
 # plot!(tmp[:,:lon],tmp[:,:lat],c=:red,w=4,leg=false)
 
 # ## 6. Interpolate Velocities
 
-dx=Γ["dx"]
+dx=Γ.dx
 uInit=[tmp[1,:lon];tmp[1,:lat]]./dx
 nSteps=Int32(tmp[end,:time]/3600)-2
 du=fill(0.0,2);
@@ -95,9 +95,9 @@ for i=1:100
 end
 
 #md plt=plot(tmpx,tmpu,label="u (interp)")
-#md plot!(Γ["XG"].f[1][1:10,1]./dx,𝑃.u0[1:10,1],marker=:o,label="u (C-grid)")
+#md plot!(Γ.XG.f[1][1:10,1]./dx,𝑃.u0[1:10,1],marker=:o,label="u (C-grid)")
 #md plot!(tmpx,tmpv,label="v (interp)")
-#md plot!(Γ["XG"].f[1][1:10,1]./dx,𝑃.v0[1:10,1],marker=:o,label="v (C-grid)")
+#md plot!(Γ.XG.f[1][1:10,1]./dx,𝑃.v0[1:10,1],marker=:o,label="v (C-grid)")
 
 # And similarly in the other direction
 
@@ -112,9 +112,9 @@ for i=1:100
 end
 
 # plt=plot(tmpx,tmpu,label="u (interp)")
-# plot!(Γ["YG"].f[1][1,1:10]./dx,𝑃.u0[1,1:10],marker=:o,label="u (C-grid)")
+# plot!(Γ.YG.f[1][1,1:10]./dx,𝑃.u0[1,1:10],marker=:o,label="u (C-grid)")
 # plot!(tmpx,tmpv,label="v (interp)")
-# plot!(Γ["YG"].f[1][1,1:10]./dx,𝑃.v0[1,1:10],marker=:o,label="v (C-grid)")
+# plot!(Γ.YG.f[1][1,1:10]./dx,𝑃.v0[1,1:10],marker=:o,label="v (C-grid)")
 
 # Compare recomputed velocities with those from `pkg/flt`
 

examples/example123.jl

@@ -131,6 +131,8 @@ function example2_setup()
    mskS=read_mds(γ.path*"hFacS",MeshArray(γ,Float32,nr))
    mskS=1.0 .+ 0.0 * mask(mskS[:,kk],NaN,0.0)
    Γ=merge(Γ,Dict("mskW" => mskW, "mskS" => mskS))
+   #convert to NamedTuple
+   Γ=(; zip(Symbol.(keys(Γ)), values(Γ))...)
 
    𝑃=FlowFields(u0[1], u1[1], v0[1], v1[1], [t0,t1])
    return 𝑃,Γ

examples/flow_fields.jl

@@ -46,17 +46,17 @@ u,v,w=solid_body_rotation(12,4)
 """
 function solid_body_rotation(np,nz)
     Γ=simple_periodic_domain(np);
-    γ=Γ["XC"].grid;
+    γ=Γ.XC.grid;
 
     #Solid-body rotation around central location ...
     i=Int(np/2+1)
-    u=-(Γ["YG"].-Γ["YG"][1][i,i])
-    v=(Γ["XG"].-Γ["XG"][1][i,i])
+    u=-(Γ.YG.-Γ.YG[1][i,i])
+    v=(Γ.XG.-Γ.XG[1][i,i])
 
     #... plus a convergent term to / from central location
     d=-0.01
-    u=u+d*(Γ["XG"].-Γ["XG"][1][i,i])
-    v=v+d*(Γ["YG"].-Γ["YG"][1][i,i])
+    u=u+d*(Γ.XG.-Γ.XG[1][i,i])
+    v=v+d*(Γ.YG.-Γ.YG[1][i,i])
 
     #Replicate u,v in vertical dimension
     uu=MeshArray(γ,γ.ioPrec,nz)
@@ -89,10 +89,10 @@ function global_ocean_circulation(;k=1,ny=2)
   p=dirname(pathof(IndividualDisplacements))
   γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
   Γ=GridLoad(γ)
-  Γ=merge(Γ,IndividualDisplacements.NeighborTileIndices_cs(Γ))
+  Γ=merge(Γ,NeighborTileIndices_cs(Γ))
 
-  func=(u -> IndividualDisplacements.update_location_llc!(u,𝐷))
-  Γ=merge(Γ,Dict("update_location!" => func))
+  func=(u -> update_location_llc!(u,𝐷))
+  Γ=merge(Γ,(; update_location! = func))
 
   #initialize u0,u1 etc
   𝑃,𝐷=set_up_FlowFields(k,Γ,ECCOclim_path);
@@ -108,19 +108,22 @@ end
 """
     OCCA_FlowFields(;backward_in_time::Bool=false)
 
-Define gridded variables and return result as Dictionary (`uvetc`).
+Define gridded variables and return result as NamedTuple
 """
 function OCCA_FlowFields(;backward_in_time::Bool=false)
 
    γ=GridSpec("PeriodicChannel",MeshArrays.GRID_LL360)
    Γ=GridLoad(γ)
-   n=length(Γ["RC"])
+   n=length(Γ.RC)
    n=5
 
-   g=Γ["XC"].grid
+   g=Γ.XC.grid
    func=(u -> IndividualDisplacements.update_location_dpdo!(u,g))
 
-   delete!.(Ref(Γ), ["hFacC", "hFacW", "hFacS","DXG","DYG","RAC","RAZ","RAS"]);
+   jj=[:hFacC, :hFacW, :hFacS, :DXG, :DYG, :RAC, :RAZ, :RAS]
+   ii=findall([!in(i,jj) for i in keys(Γ)])
+   Γ=(; zip(Symbol.(keys(Γ)[ii]), values(Γ)[ii])...)
+
    backward_in_time ? s=-1.0 : s=1.0
    s=Float32(s)
 
@@ -153,8 +156,8 @@ function OCCA_FlowFields(;backward_in_time::Bool=false)
 #   𝑆=read(rd(fileIn,"salt",n),MeshArray(γ,Float64,n))
 
    for i in eachindex(u)
-      u[i]=u[i]./Γ["DXC"][1]
-      v[i]=v[i]./Γ["DYC"][1]
+      u[i]=u[i]./Γ.DXC[1]
+      v[i]=v[i]./Γ.DYC[1]
    end
 
    for i in eachindex(u)
@@ -165,18 +168,18 @@ function OCCA_FlowFields(;backward_in_time::Bool=false)
    end
 
    for i in eachindex(w)
-      w[i]=w[i]./Γ["DRC"][min(i[2]+1,n)]
+      w[i]=w[i]./Γ.DRC[min(i[2]+1,n)]
       w[i]=circshift(w[i],[-180 0])
    end
 
-   tmpx=circshift(Γ["XC"][1],[-180 0])
+   tmpx=circshift(Γ.XC[1],[-180 0])
    tmpx[1:180,:]=tmpx[1:180,:] .- 360.0
-   Γ["XC"][1]=tmpx
+   Γ.XC[1]=tmpx
 
-   tmpx=circshift(Γ["XG"][1],[-180 0])
+   tmpx=circshift(Γ.XG[1],[-180 0])
    tmpx[1:180,:]=tmpx[1:180,:] .- 360.0
-   Γ["XG"][1]=tmpx
-   Γ["Depth"][1]=circshift(Γ["Depth"][1],[-180 0])
+   Γ.XG[1]=tmpx
+   Γ.Depth[1]=circshift(Γ.Depth[1],[-180 0])
 
    t0=0.0; t1=86400*366*2.0;
 
@@ -192,8 +195,8 @@ function OCCA_FlowFields(;backward_in_time::Bool=false)
 
    𝑃=FlowFields(u,u,v,v,w,w,[t0,t1],func)
 
-   𝐷 = (θ0=θ, θ1=θ, XC=exchange(Γ["XC"]), YC=exchange(Γ["YC"]), 
-   RF=Γ["RF"], RC=Γ["RC"],ioSize=(360,160,n))
+   𝐷 = (θ0=θ, θ1=θ, XC=exchange(Γ.XC), YC=exchange(Γ.YC), 
+   RF=Γ.RF, RC=Γ.RC,ioSize=(360,160,n))
 
    return 𝑃,𝐷,Γ
 
@@ -254,7 +257,6 @@ end
 function test2_periodic_domain(np = 12, nq = 12)
     #domain and time parameters
     Γ = simple_periodic_domain(np, nq)
-    Γ = IndividualDisplacements.dict_to_nt(Γ)
 
     u = 0.1 ./ Γ.DXC
     v = 0.3 ./ Γ.DYC

examples/helper_functions.jl

@@ -39,7 +39,7 @@ end
     isosurface(θ,T,z)
 
 ```
-isosurface(𝐼.𝑃.θ0,15,Γ["RC"])
+isosurface(𝐼.𝑃.θ0,15,Γ.RC)
 ```    
 """
 function isosurface(θ,T,z)
@@ -49,17 +49,17 @@ function isosurface(θ,T,z)
         for k=1:nr-1
             i=findall(isnan.(d[j]).&(θ[j,k].>T).&(θ[j,k+1].<=T))
             a=(θ[j,k][i] .- T)./(θ[j,k][i] .- θ[j,k+1][i])
-            d[j][i]=(1 .- a).*Γ["RC"][k] + a.*Γ["RC"][k+1]
+            d[j][i]=(1 .- a).*Γ.RC[k] + a.*Γ.RC[k+1]
             i=findall(isnan.(d[j]).&(θ[j,k].<=T).&(θ[j,k+1].>T))
             a=(θ[j,k+1][i] .- T)./(θ[j,k+1][i] .- θ[j,k][i])
-            d[j][i]=(1 .- a).*Γ["RC"][k+1] + a.*Γ["RC"][k]
+            d[j][i]=(1 .- a).*Γ.RC[k+1] + a.*Γ.RC[k]
         end
     end
     return d
 end
 
 """
-    set_up_𝑃(k::Int,t::Float64,Γ::Dict,pth::String)
+    set_up_𝑃(k::Int,t::Float64,Γ::NamedTuple,pth::String)
 
 Define `FlowFields` data structure (𝑃) along with ancillary variables (𝐷)
 for the specified grid (`Γ` dictionnary), vertical level (`k`), and 
@@ -68,24 +68,24 @@ file location (`pth`).
 _Note: the initial implementation approximates month durations to 
 365 days / 12 months for simplicity and sets 𝑃.𝑇 to [-mon/2,mon/2]_
 """
-function set_up_FlowFields(k::Int,Γ::Dict,pth::String)
-    XC=exchange(Γ["XC"]) #add 1 lon point at each edge
-    YC=exchange(Γ["YC"]) #add 1 lat point at each edge
-    iDXC=1. ./Γ["DXC"]
-    iDYC=1. ./Γ["DYC"]
-    γ=Γ["XC"].grid
+function set_up_FlowFields(k::Int,Γ::NamedTuple,pth::String)
+    XC=exchange(Γ.XC) #add 1 lon point at each edge
+    YC=exchange(Γ.YC) #add 1 lat point at each edge
+    iDXC=1. ./Γ.DXC
+    iDYC=1. ./Γ.DYC
+    γ=Γ.XC.grid
     mon=86400.0*365.0/12.0
-    func=Γ["update_location!"]
+    func=Γ.update_location!
 
     if k==0
-        msk=Γ["hFacC"]
+        msk=Γ.hFacC
         (_,nr)=size(msk)
         𝑃=FlowFields(MeshArray(γ,Float64,nr),MeshArray(γ,Float64,nr),
         MeshArray(γ,Float64,nr),MeshArray(γ,Float64,nr),
         MeshArray(γ,Float64,nr+1),MeshArray(γ,Float64,nr+1),
         [-mon/2,mon/2],func)
     else
-        msk=Γ["hFacC"][:, k]
+        msk=Γ.hFacC[:, k]
         𝑃=FlowFields(MeshArray(γ,Float64),MeshArray(γ,Float64),
         MeshArray(γ,Float64),MeshArray(γ,Float64),[-mon/2,mon/2],func)    
     end
@@ -94,8 +94,7 @@ function set_up_FlowFields(k::Int,Γ::Dict,pth::String)
          XC=XC, YC=YC, iDXC=iDXC, iDYC=iDYC,
          k=k, msk=msk, θ0=similar(msk), θ1=similar(msk))
 
-    tmp = IndividualDisplacements.dict_to_nt(IndividualDisplacements.NeighborTileIndices_cs(Γ))
-    𝐷 = merge(𝐷 , tmp)
+    𝐷 = merge(𝐷 , NeighborTileIndices_cs(Γ))
 
     return 𝑃,𝐷
 end
@@ -166,7 +165,7 @@ function update_FlowFields!(𝑃::𝐹_MeshArray3D,𝐷::NamedTuple,t::Float64)
     m1=mod(m1,12)
     m1==0 ? m1=12 : nothing
 
-    (_,nr)=size(𝐷.Γ["hFacC"])
+    (_,nr)=size(𝐷.Γ.hFacC)
 
     (U,V)=read_velocities(𝑃.u0.grid,m0,𝐷.pth)
     u0=U; v0=V
@@ -198,9 +197,9 @@ function update_FlowFields!(𝑃::𝐹_MeshArray3D,𝐷::NamedTuple,t::Float64)
     𝑃.v1[:,:]=v1[:,:]
     for k=1:nr
         tmpw=exchange(-w0[:,k],1)
-        𝑃.w0[:,k]=tmpw./𝐷.Γ["DRC"][k]
+        𝑃.w0[:,k]=tmpw./𝐷.Γ.DRC[k]
         tmpw=exchange(-w1[:,k],1)
-        𝑃.w1[:,k]=tmpw./𝐷.Γ["DRC"][k]
+        𝑃.w1[:,k]=tmpw./𝐷.Γ.DRC[k]
     end
     𝑃.w0[:,1]=0*exchange(-w0[:,1],1)
     𝑃.w1[:,1]=0*exchange(-w1[:,1],1)

examples/recipes_plots.jl

@@ -108,7 +108,7 @@ Compute Ocean depth logarithm on regular grid.
 function OceanDepthLog(Γ)
     lon=[i for i=20.:2.0:380., j=-79.:2.0:89.]
     lat=[j for i=20.:2.0:380., j=-79.:2.0:89.]
-    DL=interp_to_lonlat(Γ["Depth"],Γ,lon,lat)
+    DL=interp_to_lonlat(Γ.Depth,Γ,lon,lat)
     DL[findall(DL.<0)].=0
     DL=transpose(log10.(DL))
     DL[findall((!isfinite).(DL))].=NaN

examples/worldwide/three_dimensional_ocean.jl

@@ -42,11 +42,11 @@ To convert from longitude,latitude here we take advantage of the regularity
 of the 1 degree grid being used -- for a more general alternative, see the 
 global ocean example.
 """
-function initial_positions(Γ::Dict, nf=10000, lon_rng=(-160.0,-159.0), lat_rng=(30.0,31.0))
+function initial_positions(Γ::NamedTuple, nf=10000, lon_rng=(-160.0,-159.0), lat_rng=(30.0,31.0))
    lon=lon_rng[1] .+(lon_rng[2]-lon_rng[1]).*rand(nf)
    lat=lat_rng[1] .+(lat_rng[2]-lat_rng[1]).*rand(nf)
-   x=lon .+ (21. - Γ["XC"][1][21,1])
-   y=lat .+ (111. - Γ["YC"][1][1,111])
+   x=lon .+ (21. - Γ.XC[1][21,1])
+   y=lat .+ (111. - Γ.YC[1][1,111])
    return x,y
 end
 

src/IndividualDisplacements.jl

@@ -13,7 +13,6 @@ include("API.jl")
 include("compute.jl")
 include("data_wrangling.jl")
 include("read.jl")
-include("update_locations.jl")
 
 export Individuals, ∫!
 export FlowFields, convert_to_FlowFields

src/data_wrangling.jl

@@ -14,7 +14,7 @@ function convert_to_FlowFields(U::Array{T,2},V::Array{T,2},t1::T) where T
     np,nq=size(U)
     Γ=simple_periodic_domain(np,nq)
 
-    g=Γ["XC"].grid
+    g=Γ.XC.grid
     u=MeshArray(g,[U])
     v=MeshArray(g,[V])
     (u,v)=exchange(u,v,1)
@@ -187,8 +187,8 @@ lat=[j for i=20.:20.0:380., j=-70.:10.0:70.]
 (f,i,j,w,_,_,_)=InterpolationFactors(𝐷.Γ,vec(lon),vec(lat))
 IntFac=(lon=lon,lat=lat,f=f,i=i,j=j,w=w)
 
-tmp1=interp_to_lonlat(𝐷.Γ["Depth"],𝐷.Γ,lon,lat)
-tmp2=interp_to_lonlat(𝐷.Γ["Depth"],IntFac)
+tmp1=interp_to_lonlat(𝐷.Γ.Depth,𝐷.Γ,lon,lat)
+tmp2=interp_to_lonlat(𝐷.Γ.Depth,IntFac)
 
 ref=[5896. 5896.]
 prod(isapprox.([maximum(tmp1) maximum(tmp2)],ref,atol=1.0))
@@ -198,7 +198,7 @@ prod(isapprox.([maximum(tmp1) maximum(tmp2)],ref,atol=1.0))
 true
 ```
 """
-function interp_to_lonlat(X::MeshArray,Γ::Dict,lon,lat)
+function interp_to_lonlat(X::MeshArray,Γ::NamedTuple,lon,lat)
     (f,i,j,w,_,_,_)=InterpolationFactors(Γ,vec(lon),vec(lat))
     return reshape(Interpolate(X,f,i,j,w),size(lon))
 end