API documentation

AbstractMeshArray{T, N}

Subtype of AbstractArray{T, N}

gcmarray{T, N, AT}

gcmarray data structure. Available constructors:

gcmarray{T,N,AT}(grid::gcmgrid,meta::varmeta,f::Array{AT,N},
         fSize::Array{NTuple{N, Int}},fIndex::Array{Int,1},v::String)

gcmarray(grid::gcmgrid,f::Array{Array{T,2},N}) where {T,N}
gcmarray(grid::gcmgrid,f::Array{Array{T,N},1}) where {T,N}

gcmarray(grid::gcmgrid,fSize::Array{NTuple{N, Int}},fIndex::Array{Int,1})
gcmarray(<same as above>,n3::Int)
gcmarray(<same as above>,n3::Int,n4::Int)

gcmarray(grid::gcmgrid)
gcmarray(grid::gcmgrid,::Type{T})
gcmarray(grid::gcmgrid,::Type{T},n3::Int)
gcmarray(grid::gcmgrid,::Type{T},n3::Int,n4::Int)

gcmarray(A::Array{T,N};meta::varmeta=defaultmeta)

gcmgrid

gcmgrid data structure. Available constructors:

gcmgrid(path::String, class::String, nFaces::Int, 
        fSize::Array{NTuple{2, Int},1},
        ioSize::Union{NTuple{2, Int},Array{Int64,2}},
        ioPrec::Type, read::Function, write::Function)

The class can be set to "LatLonCap", "CubeSphere", "PeriodicChannel", "PeriodicDomain".

For example, A periodic channel (periodic in the x direction) of size 360 by 160, can be defined as follows.

pth=MeshArrays.GRID_LL360
class="PeriodicChannel"
ioSize=(360, 160)
ioPrec=Float32

γ=gcmgrid(pth, class, 1, [ioSize], ioSize, ioPrec, read, write)

Γ=GridLoad(γ)    

Please refer to GridSpec and UnitGrid for more info related to gcmgrid options.

varmeta

varmeta data structure. By default, unit is missing (non-dimensional), position is fill(0.5,3) (cell center), time is missing, and name / long_name is unknown.

Available constructors:

varmeta(unit::Union{Unitful.Units,Number,Missing},position::Array{Float64,1},
        time::Union{DateTime,Missing,Array{DateTime,1}},name::String,long_name::String)

And:

defaultmeta = varmeta(missing,fill(0.5,3),missing,"unknown","unknown")

gridpath

gridpath data structure.

using MeshArrays
γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
Γ=GridLoad(γ;option="light")
lons=[-68 -63]; lats=[-54 -66]; name="Drake Passage"
Trsct=Transect(name,lons,lats,Γ)

gridmask

gridmask data structure.

G,M,files=Integration.example()
M

gcmvector{T, N}

gcmvector data structure that can be used for subsetting and indexing into a gcmarray.

gcmvector{T,N}(grid::gcmgrid,f::Array{Array{T,1},N},
         fSize::Array{NTuple{N, Int}},fIndex::Array{Int,1})

gcmfaces{T, N}

gcmfaces data structure. Available constructors:

gcmfaces{T,N}(grid::gcmgrid,f::Array{Array{T,N},1},
         fSize::Array{NTuple{N, Int}}, aSize::NTuple{N,Int})

gcmfaces(grid::gcmgrid,v1::Array{Array{T,N},1}) where {T,N}
gcmfaces(grid::gcmgrid,::Type{T},
         fSize::Array{NTuple{N, Int}}, aSize::NTuple{N,Int}) where {T,N}

gcmfaces(grid::gcmgrid)
gcmfaces(grid::gcmgrid,::Type{T})
gcmfaces(grid::gcmgrid,::Type{T},n3::Int)

GridSpec(category="PeriodicDomain",path=tempname(); np=nothing, ID=:unknown)

• Select one of the pre-defined grids either by ID (keyword) or by category.
• Return the corresponding gcmgrid specification, including the path where grid files can be accessed (path).

1. selection by ID

• :LLC90
• :CS32
• :LLC270
• :onedegree
• :default

Example:

using MeshArrays
g = GridSpec(ID=:LLC90)

note : the path to these fully supported grids are handled internally in MeshArrays.jl.

1. by category and path

• "PeriodicDomain"
• "PeriodicChannel"
• "CubeSphere"
• "LatLonCap"`

Examples:

using MeshArrays
g = GridSpec()
g = GridSpec("PeriodicChannel",MeshArrays.GRID_LL360)
g = GridSpec("CubeSphere",MeshArrays.GRID_CS32)
g = GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
isa(g,gcmgrid)

1. by category and path with the np argument

• np is the number of grid points in x or y for the LatLonCap and CubeSphere tiles.

np defaults to 90 for LatLonCap and 32 for CubeSphere, and so must be included to access the LLC270 grid with the category argument.

Examples:

using MeshArrays
g = GridSpec("LatLonCap",MeshArrays.GRID_LLC90,np=90)
g = GridSpec("LatLonCap",MeshArrays.GRID_LLC270,np=270)
g = GridSpec("CubeSphere",MeshArrays.GRID_CS32,np=32)
isa(g,gcmgrid)

GridLoad(γ=GridSpec(); ID=:default, option=:minimal)

• Return a NamedTuple of grid variables read from files located in γ.path (see ?GridSpec).
• option :
  • option=:minimal (default) to get only grid cell center positions (XC, YC).
  • option=:light to get a complete set of 2D grid variables.
  • option=:full to get a complete set of 2D & 3D grid variables.

Based on the MITgcm naming convention, grid variables are:

• XC, XG, YC, YG, AngleCS, AngleSN, hFacC, hFacS, hFacW, Depth.
• RAC, RAW, RAS, RAZ, DXC, DXG, DYC, DYG.
• DRC, DRF, RC, RF (one-dimensional)

MITgcm documentation :

https://mitgcm.readthedocs.io/en/latest/algorithm/algorithm.html#spatial-discretization-of-the-dynamical-equations

using MeshArrays
γ = GridSpec(ID=:LLC90)
Γ = GridLoad(γ;option="full")

isa(Γ.XC,MeshArray)

GridLoadVar(nam::String,γ::gcmgrid)

Return a grid variable read from files located in γ.path (see ?GridSpec, ?GridLoad).

Based on the MITgcm naming convention, grid variables are:

• XC, XG, YC, YG, AngleCS, AngleSN, hFacC, hFacS, hFacW, Depth.
• RAC, RAW, RAS, RAZ, DXC, DXG, DYC, DYG.
• DRC, DRF, RC, RF (one-dimensional)

MITgcm documentation :

https://mitgcm.readthedocs.io/en/latest/algorithm/algorithm.html#spatial-discretization-of-the-dynamical-equations

using MeshArrays

γ = GridSpec("CubeSphere",MeshArrays.GRID_CS32)
XC = GridLoadVar("XC",γ)

isa(XC,MeshArray)

periodic_domain(np::Integer,nq=missing)

Set up a simple periodic domain of size np x nq

using MeshArrays
np=16 #domain size is np x np
Γ=Grids_simple.periodic_domain(np)
isa(Γ.XC,MeshArray)

grid_factors(xy::NamedTuple)

xy=Grids_simple.xy_OISST()
gr=Grids_simple.grid_factors(xy)

grid_add_z(G,depth,mask)

xy=Grids_simple.xy_IAP()
gr=Grids_simple.grid_factors(xy)

dep=[10 100 1000]; msk=ones(gr[:XC].fSize[1]...,3)
gr=Grids_simple.grid_add_z(gr,dep,msk)

UnitGrid(γ::gcmgrid;option="minimal")

Generate a unit grid, where every grid spacing and area is 1, according to γ.

UnitGrid(ioSize, tileSize; option="minimal")

Generate a unit grid, where every grid spacing and area is 1, according to ioSize, tileSize.

Since ioSize, tileSize defines a one to one mapping from global to tiled array, here we use read_tiles, write_tiles instead of the default read, write. And we overwrite local XC,YC etc accordingly with global XC,YC etc .

Tiles(γ::gcmgrid,ni::Int,nj::Int)

Define sudomain tiles of size ni,nj. Each tile is defined by a Dict where tile,face,i,j correspond to tile ID, face ID, index ranges.

using MeshArrays
γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
τ=Tiles(γ,30,30)

isa(τ[1],NamedTuple)

Tiles(τ::Array{Dict},x::MeshArray)

Return an Array of tiles which cover x according to tile partition τ.

using MeshArrays
γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
d=γ.read(γ.path*"Depth.data",MeshArray(γ,γ.ioPrec))
τ=Tiles(γ,30,30)
td=Tiles(τ,d)

D=similar(d)
Tiles!(τ,td,D)

isa(td[1],Array)

Tiles!(τ::Array,tx::Array,x::MeshArrays)

Map tiles in tx according to tile partition τ into x.

exchange(fld::MeshArray)

Exchange / transfer data between neighboring arrays. Other methods are

exchange(fld::MeshArray,N::Integer)
exchange(u::MeshArray,v::MeshArray)
exchange(u::MeshArray,v::MeshArray,N::Integer)

read(fil::String,x::MeshArray)

Read array from file and return as a MeshArray.

The second argument (MeshArray or gcmgrid) provides the grid specifications (x.grid.ioSize). ```

read(xx::Array,γ::gcmgrid)

Reformat Array data into a MeshArray shaped after γ.

read(xx::Array,x::MeshArray)

Reformat Array data into a MeshArray similar to x.

read!(xx::Array,x::MeshArray)

Reformat array into MeshArray and write into x.

write(fil::String,x::MeshArray)

Write MeshArray to binary file. Other methods:

write(xx::Array,x::MeshArray) #to Array

interpolation_setup(fil::String)

Read e.g. interp_coeffs_halfdeg.jld2

fil=joinpath(tempdir(),"interp_coeffs_halfdeg.jld2")
λ=MeshArrays.interpolation_setup(fil)

interpolation_setup(;Γ,lon,lat,filename)

Download or recompute interpolation coefficients.

• λ=interpolation_setup() to download "interpcoeffshalfdeg.jld2"
• λ=interpolation_setup(Γ=Γ) to recompute interpolation to lon,lat
• λ=interpolation_setup(Γ=Γ,lon=lon,lat=lat,filename=filename) to recompute interpolation to lon,lat and save to location and file filename, defaulting to tempname()*"_interp_coeffs.jld2".

Interpolate(z_in::MeshArray,f,i,j,w)

using MeshArrays

γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
Γ=GridLoad(γ; option="full")

lon=[i for i=-179.5:1.0:179.5, j=-89.5:1.0:89.5]
lat=[j for i=-179.5:1.0:179.5, j=-89.5:1.0:89.5]
(f,i,j,w,j_f,j_x,j_y)=InterpolationFactors(Γ,vec(lon),vec(lat))
DD=Interpolate(Γ.Depth,f,i,j,w)

using CairoMakie
heatmap(vec(lon[:,1]),vec(lat[1,:]),DD,colorrange=(0.,6000.))

Interpolate(z_in::MeshArray,f,i,j,w)

using MeshArrays

γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
Γ=GridLoad(γ; option="full")

lon=[i for i=-179.5:1.0:179.5, j=-89.5:1.0:89.5]
lat=[j for i=-179.5:1.0:179.5, j=-89.5:1.0:89.5]
(f,i,j,w,j_f,j_x,j_y)=InterpolationFactors(Γ,vec(lon),vec(lat))
L=(lon=lon, lat=lat, f=f, i=i, j=j, w=w)

using CairoMakie
heatmap(Interpolate(Γ.Depth,L)...,colorrange=(0.,6000.))

or

heatmap(Γ.Depth,interpolation=L,colorrange=(0.,6000.))

InterpolationFactors(Γ,lon::Array{T,1},lat::Array{T,1})

Compute interpolation coefficients etc from grid Γ to lon,lat

using MeshArrays
γ=GridSpec("CubeSphere",MeshArrays.GRID_CS32)
Γ=GridLoad(γ; option="full")
lon=collect(45.:0.1:46.); lat=collect(60.:0.1:61.)
(f,i,j,w,j_f,j_x,j_y)=InterpolationFactors(Γ,lon,lat)
YC=Interpolate(Γ.YC,f,i,j,w)
extrema(i)==(9,10)

StereographicProjection(XC0,YC0,XC,YC)

Apply stereographic projection that puts XC0,YC0 at 0.0,0.0 to target point(s) XC,YC

lon=collect(45.:0.1:46.); lat=collect(60.:0.1:61.)
x,y=StereographicProjection(45.,60.,lon,lat)

knn(xgrid,ygrid::MeshArray,x,y::Array{T,1},k::Int)

Find k nearest neighbors to each point in x,y on xgrid,ygrid

lon=collect(0.1:0.5:2.1); lat=collect(0.1:0.5:2.1);
(f,i,j,c)=knn(Γ.XC,Γ.YC,lon,lat)

curl(u::MeshArray,v::MeshArray,Γ::NamedTuple)

Compute curl of a velocity field.

convergence(uFLD::MeshArray,vFLD::MeshArray)

Compute convergence of a vector field

gradient(inFLD::MeshArray,Γ::NamedTuple)

Compute spatial derivatives. Other methods:

gradient(inFLD::MeshArray,Γ::NamedTuple,doDIV::Bool)
gradient(inFLD::MeshArray,iDXC::MeshArray,iDYC::MeshArray)

ScalarPotential(TrspCon)

Scalar potential inversion.

TrspPot=ScalarPotential(TrspCon)

VectorPotential(TrspX,TrspY,Γ,method::Int=1)

Vector potential inversion.

TrspPot=ScalarPotential(TrspCon)

ThroughFlow(VectorField,IntegralPath,Γ::NamedTuple)

Compute transport through an integration path

UVtoTransport(U,V,G::NamedTuple)

Convert e.g. velocity (m/s) to transport (m^3/s) by multiplying by DRF and DXG,DYG.

UVtoUEVN(u,v,G::NamedTuple)

1. Interpolate to grid cell centers (uC,vC)
2. Convert to Eastward/Northward components (uE,vN)

Note: land masking u,v with NaNs preemptively can be adequate.

loops(mask::gridmask; files=String[], var=:THETA, rd=read)

begin
  @everywhere using MeshArrays, MITgcm
  @everywhere rd(F,var,tim,tmp)=read(read_mdsio(F,var),tmp)
  @everywhere G,M,files=Integration.example()
    #,regions=(30,10),depths=Integration.DEPTHS)
end;

H=Integration.loops(M,files=files,rd=rd)
# Hbis=Integration.streamlined_loop(M,files=files,rd=rd)

and to save results:

using JLD2; output_path="test.jld2"
jldsave(output_path; depths=M.depths, integral=H, volume=vol, name=M.names)

where vol is calculated as follows:

#option=:loops
allones=1.0 .+0*G.hFacC
M.tmp2d.=M.v_int[1](allones)
vol=[b(M.tmp2d) for b in M.h_sum]

or

#option=:streamlined_loop
allones=1.0 .+0*G.hFacC
vol=[b(allones) for b in M.h_sum]

LatitudeCircles(LatValues,Γ::NamedTuple; format=:gridpath, range=(0.0,360.0))

Compute integration paths that follow latitude circles, within the specified longitude range.

Transect(name,lons,lats,Γ; segment=:short, format=:gridpath)

Compute integration paths that follow a great circle between two geolocations give by lons, lats.

ocean_basins()

Mask of ocean basins on ECCO grid (LLC90).

basins=demo.ocean_basins()
AtlExt=extended_basin(basins,:Atl)

extended_basin(basins,ID=:Atl)

Consolidate basins mask to include marginal seas.

note : has only be tested on the ECCO grid (LLC90).

basins=demo.ocean_basins()
AtlExt=extended_basin(basins,:Atl)

isosurface(θ,T,Γ)

Depth of isosurface θ=T

mydatadep(nam="countries_shp1")

Download data dependency with predefined name; currently :

"countries_shp1"
"countries_geojson1"
"basemap_jpg1"
"interp_halfdeg"
"oceans_geojson1"