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.

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

γ=gcmgrid(path, 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.

path=MeshArrays.Dataset("GRID_LLC90")
γ=GridSpec("LatLonCap",path)
Γ=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="default",
    path=tempname(); np=nothing, ID=:unknown)

• Select one of the pre-defined grids
  • either by category parameter (e.g. "default" or "ones")
  • or by ID keyword via MeshArrays.GridSpec_default or MeshArrays.GridSpec_MITgcm)
• Return the corresponding gcmgrid
  • including in path either a path to grid files or a placeholder like _ones or _default.

using MeshArrays
γ=GridSpec()

GridSpec_default(xy=NamedTuple(), nFaces=1; ID=:unknown)

• Select one of the pre-defined grids
  • or by providing xy parameter (a NamedTuple that includes :xc, :yc, :xg, :yg)
  • either by ID keyword (e.g., :OISST, :Oscar, or :IAP by default)
• Return the corresponding gcmgrid
  • incl. path="_default"

Example:

using MeshArrays
a = MeshArrays.GridSpec_default(ID=:OISST)
b = MeshArrays.GridSpec_default(Grids_simple.xy_OISST())

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

• Select one of the pre-defined grids
  • or by category parameter
  • either by ID keyword
• Return the corresponding gcmgrid
  • incl. path to grid files (path).

Input files for these fully supported grids get downloaded internally via MeshArrays.Dataset if needed.

1. selection by ID

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

Example:

using MeshArrays
g = GridSpec(ID=:LLC90)

1. by category and path

• "PeriodicDomain" (default)
• "PeriodicChannel"
• "CubeSphere"
• "LatLonCap"`

Examples:

using MeshArrays
g = GridSpec("CubeSphere",MeshArrays.Dataset("GRID_CS32"))
isa(g,gcmgrid)

When the np keyword argument is specified then it is used to set the number of points in each grid subdomain.

Example:

using MeshArrays
g = GridSpec("LatLonCap",np=270)
isa(g,gcmgrid)

GridSpec_ones(grTp,nF,nP)

Define grid where all variables will be set to one via GridLoad_ones.

using MeshArrays
γ=MeshArrays.GridSpec_ones("CubeSphere",6,20);
Γ=MeshArrays.GridLoad_ones(γ;option="full")

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

• Return a NamedTuple of grid variables read from files located in γ.path (see ?GridSpec).
• if ID is specified then call GridSpec(ID=ID) to override γ parameter.
• 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.

For grid variables, we follow the MITgcm naming convention. Grid variables thus typically include :

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

For additional detail please refer to the 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)

GridLoad_default(γ=GridSpec())

MeshArrays.GridSpec_default(ID=:IAP)

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

Generate a unit grid, where every grid spacing and area is 1, as specified by γ obtained from GridSpec_ones.

using MeshArrays
γ=MeshArrays.GridSpec_ones("CubeSphere",6,20);
Γ=MeshArrays.GridLoad_ones(γ;option="full")

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.Dataset("GRID_CS32"))
XC = GridLoadVar("XC",γ)

isa(XC,MeshArray)

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.Dataset("GRID_LLC90"))
τ=Tiles(γ,30,30)

isa(τ[1],NamedTuple)

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

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

using MeshArrays
γ=GridSpec("LatLonCap",MeshArrays.Dataset("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::AbstractMeshArray)

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

exchange(fld::AbstractMeshArray)

Exchange / transfer data between neighboring arrays. Other methods are

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

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

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::AbstractMeshArray)

Reformat Array data into a MeshArray similar to x.

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

Reformat array into MeshArray and write into x.

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

Write MeshArray to binary file. Other methods:

write(xx::Array,x::AbstractMeshArray) #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::AbstractMeshArray,f,i,j,w)

using MeshArrays
path=MeshArrays.Dataset("GRID_LLC90")
γ=GridSpec("LatLonCap",path)
Γ=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::AbstractMeshArray,f,i,j,w)

using MeshArrays
path=MeshArrays.Dataset("GRID_LLC90")
γ=GridSpec("LatLonCap",path)
Γ=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
path=MeshArrays.Dataset("GRID_CS32")
γ=GridSpec("CubeSphere",path)
Γ=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::AbstractMeshArray,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::AbstractMeshArray,v::AbstractMeshArray,Γ::NamedTuple)

Compute curl of a velocity field.

convergence(uFLD::AbstractMeshArray,vFLD::AbstractMeshArray)

Compute convergence of a vector field

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

Compute spatial derivatives. Other methods:

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

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]

Dataset(d::String; do_read=true, verbose=false)

• get folder name for dataset d.
• download dataset if needed.
• read content if do_read

mydatadep(nam="countries_shp1")

Download data dependency with predefined name; currently :

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

struct NamedPolygon geometry::Vector{Tuple{Float32, Float32}} name::String points::Vector{} end

struct polyarray name::String f::Array{NamedPolygon} end

using MeshArrays, GeoJSON, DataDeps
pol=MeshArrays.Dataset("oceans_geojson1")

write(pol,tempname()*".json")

using CairoMakie
lines(pol)

to_Polygon(pa::polyarray)

Convert polyarray into a GI.Polygon array.

pol_P,nams=MeshArrays.to_Polygon(pol)

 to_polyarray(pol)

Convert polygon data into a polyarray.

fil="countries.geojson"
MeshArrays.to_polyarray(GeoJSON.read(fil))

within_pol(pol; ID=1)

Generate a name,rule pair to test if location lon,lat is within polygon pol[ID].geometry.

using MeshArrays, DataDeps, GeoJSON, GeometryOps
fil=MeshArrays.Dataset("oceans_geojson1",do_read=false)
pol=MeshArrays.read_json(fil,format=:GeoJSON)

name,rule=MeshArrays.within_pol(pol; ID=11)
rule(-30,40)

read_json(fil; format=:polyarray)

Call GeoJSON.read and return polyarray (default)

• format :polyarray(default): returnpolyarray`
• format :coord : convert to GI.coordinates.
• format :GeoJSON : return FeatureCollection (vector of name,geom pairs)

import MeshArrays, DataDeps, GeoJSON
pol=MeshArrays.Dataset("oceans_geojson1")

read_shp(fil; format=:polyarray)

Call Shapefile.Table and Shapefile.shapes and return polyarray (default).

• format :polyarray (default) : return polyarray
• format :Shapefile : return vector of name,geometry` named tuples.
• format :coord : convert to GI.coordinates.

using MeshArrays, DataDeps, Shapefile
fil=MeshArrays.Dataset("countries_shp1")

using CairoMakie
lines(pol)

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