16A_run_divand_with_cv_jl.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the DIVAnd analysis, using DIVAnd cross-validation to estimate the correlation length and epsilon2 input parameters. Use the domain of the observations as the input domain (not the domain of the mask)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using DIVAnd\n",
    "using CSV\n",
    "using NCDatasets\n",
    "using DataFrames\n",
    "using Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/home/hourstonh/Documents/climatology/data/value_vs_depth/16_diva_analysis/pmn/\""
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Choose depth/year/season and set path names\n",
    "\n",
    "var_name = \"Oxy\"\n",
    "standard_depth = 0\n",
    "year = 1999\n",
    "szn = \"OND\"\n",
    "# subsamp_interval = 1\n",
    "\n",
    "# obs_dir = string(\"C:\\\\Users\\\\HourstonH\\\\Documents\\\\NEP_climatology\\\\data\\\\\",\n",
    "#     \"value_vs_depth\\\\14_sep_by_sl_and_year\\\\\")\n",
    "\n",
    "# mask_dir = string(\"C:\\\\Users\\\\HourstonH\\\\Documents\\\\NEP_climatology\\\\data\\\\\",\n",
    "#     \"value_vs_depth\\\\16_diva_analysis\\\\masks\\\\\")\n",
    "\n",
    "# pmn_dir = string(\"C:\\\\Users\\\\HourstonH\\\\Documents\\\\NEP_climatology\\\\data\\\\\",\n",
    "#     \"value_vs_depth\\\\16_diva_analysis\\\\pmn\\\\\")\n",
    "\n",
    "# output_dir = string(\"C:\\\\Users\\\\HourstonH\\\\Documents\\\\NEP_climatology\\\\data\\\\\",\n",
    "#     \"value_vs_depth\\\\16_diva_analysis\\\\\")\n",
    "\n",
    "# Paths for Linux:\n",
    "obs_dir = string(\"/home/hourstonh/Documents/climatology/data/\",\n",
    "    \"value_vs_depth/14_sep_by_sl_and_year/\")\n",
    "\n",
    "mask_dir = string(\"/home/hourstonh/Documents/climatology/data/\",\n",
    "    \"value_vs_depth/16_diva_analysis/masks/\")\n",
    "\n",
    "pmn_dir = string(\"/home/hourstonh/Documents/climatology/data/\",\n",
    "    \"value_vs_depth/16_diva_analysis/pmn/\")\n",
    "\n",
    "# years = [1991:1:2020;]\n",
    "# println(years)\n",
    "# szns = [\"JFM\", \"AMJ\", \"JAS\", \"OND\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/hourstonh/Documents/climatology/data/value_vs_depth/14_sep_by_sl_and_year/Oxy_0m_1999_OND.csv\n",
      "0.7018335999999991 15.941333599999998\n"
     ]
    }
   ],
   "source": [
    "# Open the required files\n",
    "\n",
    "# Read in standard level data file\n",
    "obs_filename = string(obs_dir, var_name, \"_\", standard_depth, \"m_\", year, \n",
    "    \"_\", szn, \".csv\")\n",
    "    \n",
    "println(obs_filename)\n",
    "    \n",
    "# Pipe operator to dataframe\n",
    "obs_df = CSV.File(obs_filename) |> DataFrame\n",
    "    \n",
    "if size(obs_df)[1] == 0\n",
    "    println(\"DataFrame empty -- skip\")\n",
    "end\n",
    "    \n",
    "xobs = obs_df[!, :Longitude]\n",
    "yobs = obs_df[!, :Latitude]\n",
    "vobs = obs_df[!, :SL_value]\n",
    "\n",
    "# Compute anomaly field\n",
    "vmean = mean(vobs)\n",
    "vanom = vobs .- vmean\n",
    "    \n",
    "# Calculate domain size based on the observations\n",
    "# Set first guesses for correlation length as 1/10 domain size\n",
    "# Should be in meters??\n",
    "domain_size_x_deg = maximum(xobs) - minimum(xobs)\n",
    "domain_size_y_deg = maximum(yobs) - minimum(xobs)\n",
    "# lenx_guess = DIVAnd.deg2m(domain_size_x_deg/10)  # /10\n",
    "# leny_guess = DIVAnd.deg2m(domain_size_y_deg/10)\n",
    "\n",
    "lenx_guess = domain_size_x_deg/10\n",
    "leny_guess = domain_size_y_deg/10\n",
    "\n",
    "println(lenx_guess, \" \", leny_guess)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "closed NetCDF NCDataset"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Read in mask\n",
    "mask_filename = string(mask_dir, var_name, \"_\", standard_depth, \"m_\", \n",
    "    year, \"_\", szn, \"_mask_6min.nc\")\n",
    "    \n",
    "mask_ds = Dataset(mask_filename)\n",
    "    \n",
    "mask = mask_ds[\"mask\"][:,:]  # [1:subsamp_interval:end, 1:subsamp_interval:end]\n",
    "mask = Bool.(mask)\n",
    "#     println(typeof(mask))\n",
    "\n",
    "# Open the pmn netCDF file\n",
    "pmn_filename = string(pmn_dir, \"divand_pmn_for_mask_6min.nc\")\n",
    "\n",
    "pmn_ds = Dataset(pmn_filename)\n",
    "\n",
    "# Get the inverse of the resolution of the grid\n",
    "pm = convert(Array{Float64}, pmn_ds[\"pm\"][:,:])\n",
    "pn = convert(Array{Float64}, pmn_ds[\"pn\"][:,:])\n",
    "\n",
    "# Get the 2d mesh grid longitude and latitude\n",
    "Lon2d, Lat2d = ndgrid(mask_ds[\"lon\"][:], mask_ds[\"lat\"][:])\n",
    "\n",
    "close(mask_ds)\n",
    "close(pmn_ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Set some more parameters\n",
    "\n",
    "signal_to_noise_ratio = 50.  # Default from Lu ODV session\n",
    "epsilon2_guess = 1/signal_to_noise_ratio  # 1.\n",
    "\n",
    "# Choose number of testing points around the current value of L (corlen)\n",
    "nl = 0\n",
    "\n",
    "# Choose number of testing points around the current value of epsilon2\n",
    "ne = 0\n",
    "\n",
    "# Choose cross-validation method\n",
    "# 1: full CV; 2: sampled CV; 3: GCV; 0: automatic choice between the three\n",
    "method = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "OutOfMemoryError()",
     "output_type": "error",
     "traceback": [
      "OutOfMemoryError()",
      "",
      "Stacktrace:",
      "  [1] _growend!",
      "    @ ./array.jl:888 [inlined]",
      "  [2] resize!",
      "    @ ./array.jl:1108 [inlined]",
      "  [3] sparse!(I::Vector{Int64}, J::Vector{Int64}, V::Vector{Float64}, m::Int64, n::Int64, combine::typeof(+), klasttouch::Vector{Int64}, csrrowptr::Vector{Int64}, csrcolval::Vector{Int64}, csrnzval::Vector{Float64}, csccolptr::Vector{Int64}, cscrowval::Vector{Int64}, cscnzval::Vector{Float64})",
      "    @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:914",
      "  [4] sparse(I::Vector{Int64}, J::Vector{Int64}, V::Vector{Float64}, m::Int64, n::Int64, combine::Function)",
      "    @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:767",
      "  [5] sparse",
      "    @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:956 [inlined]",
      "  [6] _sparse_wsum(sz1::Tuple{Int64, Int64}, m::Int64, v1::Float64, v2::Float64, cyclic::Bool)",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/sparse_diff.jl:34",
      "  [7] sparse_diff",
      "    @ ~/.julia/packages/DIVAnd/eDBbE/src/sparse_diff.jl:48 [inlined]",
      "  [8] oper_diff",
      "    @ ~/.julia/packages/DIVAnd/eDBbE/src/function_operator.jl:14 [inlined]",
      "  [9] DIVAnd_laplacian(operatortype::Type, mask::BitMatrix, pmn::Tuple{Matrix{Float64}, Matrix{Float64}}, nu::Tuple{Matrix{Float64}, Matrix{Float64}}, iscyclic::BitVector; coeff_laplacian::Vector{Float64})",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAnd_laplacian.jl:97",
      " [10] DIVAnd_operators(operatortype::Type, mask::BitMatrix, pmn::Tuple{Matrix{Float64}, Matrix{Float64}}, nu::Tuple{Matrix{Float64}, Matrix{Float64}}, iscyclic::BitVector, mapindex::Vector{Any}, Labs::Tuple{Matrix{Float64}, Matrix{Float64}}; coeff_laplacian::Vector{Float64})",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAnd_operators.jl:59",
      " [11] DIVAnd_background(operatortype::Type, mask::BitMatrix, pmn::Tuple{Matrix{Float64}, Matrix{Float64}}, Labs::Tuple{Matrix{Float64}, Matrix{Float64}}, alpha::Vector{Any}, moddim::Vector{Any}, scale_len::Bool, mapindex::Vector{Any}; btrunc::Vector{Any}, coeff_laplacian::Vector{Float64}, coeff_derivative2::Vector{Float64}, mean_Labs::Nothing)",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAnd_background.jl:92",
      " [12] DIVAndrun(operatortype::Type, mask::BitMatrix, pmnin::Tuple{Matrix{Float64}, Matrix{Float64}}, xiin::Tuple{Matrix{Float64}, Matrix{Float64}}, x::Tuple{Vector{Float64}, Vector{Float64}}, f::Vector{Float64}, lin::Tuple{Float64, Float64}, epsilon2::Float64; velocity::Tuple{}, primal::Bool, factorize::Bool, tol::Float64, maxit::Int64, minit::Int64, constraints::Tuple{}, inversion::Symbol, moddim::Vector{Any}, fracindex::Matrix{Float64}, alpha::Vector{Any}, keepLanczosVectors::Int64, compPC::typeof(DIVAnd_pc_none), progress::DIVAnd.var\"#275#277\", fi0::Matrix{Float64}, f0::Vector{Float64}, alphabc::Float64, scale_len::Bool, btrunc::Vector{Any}, MEMTOFIT::Float64, topographyforfluxes::Tuple{}, fluxes::Tuple{}, epsfluxes::Int64, epsilon2forfractions::Int64, RTIMESONESCALES::Tuple{}, QCMETHOD::Tuple{}, coeff_laplacian::Vector{Float64}, coeff_derivative2::Vector{Float64}, mean_Labs::Nothing)",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAndrun.jl:53",
      " [13] DIVAndrun(operatortype::Type, mask::BitMatrix, pmnin::Tuple{Matrix{Float64}, Matrix{Float64}}, xiin::Tuple{Matrix{Float64}, Matrix{Float64}}, x::Tuple{Vector{Float64}, Vector{Float64}}, f::Vector{Float64}, lin::Tuple{Float64, Float64}, epsilon2::Float64)",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAndrun.jl:45",
      " [14] #DIVAndrun#280",
      "    @ ~/.julia/packages/DIVAnd/eDBbE/src/DIVAndrun.jl:298 [inlined]",
      " [15] DIVAndrun",
      "    @ ~/.julia/packages/DIVAnd/eDBbE/src/DIVAndrun.jl:298 [inlined]",
      " [16] DIVAnd_cv(mask::BitMatrix, pmn::Tuple{Matrix{Float64}, Matrix{Float64}}, xi::Tuple{Matrix{Float64}, Matrix{Float64}}, x::Tuple{Vector{Float64}, Vector{Float64}}, f::Vector{Float64}, len::Tuple{Float64, Float64}, epsilon2::Float64, nl::Int64, ne::Int64, method::Int64; rng::Random._GLOBAL_RNG, otherargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAnd_cv.jl:136",
      " [17] DIVAnd_cv(mask::BitMatrix, pmn::Tuple{Matrix{Float64}, Matrix{Float64}}, xi::Tuple{Matrix{Float64}, Matrix{Float64}}, x::Tuple{Vector{Float64}, Vector{Float64}}, f::Vector{Float64}, len::Tuple{Float64, Float64}, epsilon2::Float64, nl::Int64, ne::Int64, method::Int64)",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAnd_cv.jl:83",
      " [18] top-level scope",
      "    @ In[6]:2",
      " [19] eval",
      "    @ ./boot.jl:360 [inlined]",
      " [20] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
      "    @ Base ./loading.jl:1116"
     ]
    }
   ],
   "source": [
    "# Run the cross-validation\n",
    "bestfactorl, bestfactore, cvval, cvvalues, x2Ddata, y2Ddata, cvinter, xi2D, yi2D = DIVAnd_cv(\n",
    "    mask, (pm, pn), (Lon2d, Lat2d), (xobs, yobs), vanom, (lenx_guess, leny_guess), \n",
    "    epsilon2_guess, nl, ne, method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply factors to guesses\n",
    "new_lenx = bestfactorl .* lenx_guess\n",
    "new_leny = bestfactorl .* leny_guess\n",
    "new_epsilon2 = bestfactore .* epsilon2_guess\n",
    "\n",
    "old_params = [lenx_guess, leny_guess, epsilon2_guess]\n",
    "new_params = [new_lenx, new_leny, new_epsilon2]\n",
    "\n",
    "for i=[1,2,3]\n",
    "    println(\"lenx:\", old_params[i], \"->\", new_params[i])\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "OutOfMemoryError()",
     "output_type": "error",
     "traceback": [
      "OutOfMemoryError()",
      "",
      "Stacktrace:",
      "  [1] _growend!",
      "    @ ./array.jl:888 [inlined]",
      "  [2] resize!",
      "    @ ./array.jl:1108 [inlined]",
      "  [3] sparse!(I::Vector{Int64}, J::Vector{Int64}, V::Vector{Float64}, m::Int64, n::Int64, combine::typeof(+), klasttouch::Vector{Int64}, csrrowptr::Vector{Int64}, csrcolval::Vector{Int64}, csrnzval::Vector{Float64}, csccolptr::Vector{Int64}, cscrowval::Vector{Int64}, cscnzval::Vector{Float64})",
      "    @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:914",
      "  [4] sparse(I::Vector{Int64}, J::Vector{Int64}, V::Vector{Float64}, m::Int64, n::Int64, combine::Function)",
      "    @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:767",
      "  [5] sparse",
      "    @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:956 [inlined]",
      "  [6] _sparse_wsum(sz1::Tuple{Int64, Int64}, m::Int64, v1::Float64, v2::Float64, cyclic::Bool)",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/sparse_diff.jl:34",
      "  [7] sparse_diff",
      "    @ ~/.julia/packages/DIVAnd/eDBbE/src/sparse_diff.jl:48 [inlined]",
      "  [8] oper_diff",
      "    @ ~/.julia/packages/DIVAnd/eDBbE/src/function_operator.jl:14 [inlined]",
      "  [9] DIVAnd_laplacian(operatortype::Type, mask::BitMatrix, pmn::Tuple{Matrix{Float64}, Matrix{Float64}}, nu::Tuple{Matrix{Float64}, Matrix{Float64}}, iscyclic::BitVector; coeff_laplacian::Vector{Float64})",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAnd_laplacian.jl:97",
      " [10] DIVAnd_operators(operatortype::Type, mask::BitMatrix, pmn::Tuple{Matrix{Float64}, Matrix{Float64}}, nu::Tuple{Matrix{Float64}, Matrix{Float64}}, iscyclic::BitVector, mapindex::Vector{Any}, Labs::Tuple{Matrix{Float64}, Matrix{Float64}}; coeff_laplacian::Vector{Float64})",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAnd_operators.jl:59",
      " [11] DIVAnd_background(operatortype::Type, mask::BitMatrix, pmn::Tuple{Matrix{Float64}, Matrix{Float64}}, Labs::Tuple{Matrix{Float64}, Matrix{Float64}}, alpha::Vector{Any}, moddim::Vector{Any}, scale_len::Bool, mapindex::Vector{Any}; btrunc::Vector{Any}, coeff_laplacian::Vector{Float64}, coeff_derivative2::Vector{Float64}, mean_Labs::Nothing)",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAnd_background.jl:92",
      " [12] DIVAndrun(operatortype::Type, mask::BitMatrix, pmnin::Tuple{Matrix{Float64}, Matrix{Float64}}, xiin::Tuple{Matrix{Float64}, Matrix{Float64}}, x::Tuple{Vector{Float64}, Vector{Float64}}, f::Vector{Float64}, lin::Tuple{Float64, Float64}, epsilon2::Float64; velocity::Tuple{}, primal::Bool, factorize::Bool, tol::Float64, maxit::Int64, minit::Int64, constraints::Tuple{}, inversion::Symbol, moddim::Vector{Any}, fracindex::Matrix{Float64}, alpha::Vector{Any}, keepLanczosVectors::Int64, compPC::typeof(DIVAnd_pc_none), progress::DIVAnd.var\"#275#277\", fi0::Matrix{Float64}, f0::Vector{Float64}, alphabc::Float64, scale_len::Bool, btrunc::Vector{Any}, MEMTOFIT::Float64, topographyforfluxes::Tuple{}, fluxes::Tuple{}, epsfluxes::Int64, epsilon2forfractions::Int64, RTIMESONESCALES::Tuple{}, QCMETHOD::Tuple{}, coeff_laplacian::Vector{Float64}, coeff_derivative2::Vector{Float64}, mean_Labs::Nothing)",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAndrun.jl:53",
      " [13] DIVAndrun(operatortype::Type, mask::BitMatrix, pmnin::Tuple{Matrix{Float64}, Matrix{Float64}}, xiin::Tuple{Matrix{Float64}, Matrix{Float64}}, x::Tuple{Vector{Float64}, Vector{Float64}}, f::Vector{Float64}, lin::Tuple{Float64, Float64}, epsilon2::Float64)",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAndrun.jl:45",
      " [14] #DIVAndrun#280",
      "    @ ~/.julia/packages/DIVAnd/eDBbE/src/DIVAndrun.jl:298 [inlined]",
      " [15] DIVAndrun(mask::BitMatrix, pmnin::Tuple{Matrix{Float64}, Matrix{Float64}}, xiin::Tuple{Matrix{Float64}, Matrix{Float64}}, x::Tuple{Vector{Float64}, Vector{Float64}}, f::Vector{Float64}, lin::Tuple{Float64, Float64}, epsilon2::Float64)",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAndrun.jl:298",
      " [16] DIVAndrunfi(::BitMatrix, ::Vararg{Any, N} where N)",
      "    @ DIVAnd ~/.julia/packages/DIVAnd/eDBbE/src/DIVAndrun.jl:157",
      " [17] top-level scope",
      "    @ In[8]:2",
      " [18] eval",
      "    @ ./boot.jl:360 [inlined]",
      " [19] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
      "    @ Base ./loading.jl:1116"
     ]
    }
   ],
   "source": [
    "# Run the analysis\n",
    "va = DIVAndrunfi(mask, (pm, pn), (Lon2d, Lat2d), (xobs, yobs), vanom,\n",
    "                 (lenx_guess, leny_guess), epsilon2_guess)\n",
    "\n",
    "# Add the output anomaly back to the mean of the observations\n",
    "vout = vmean .+ va"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export the pm, pn as a netCDF file in the same dims as the mask\n",
    "# pmn_filename = string(pmn_dir, var_name, \"_\", standard_depth, \"m_\",\n",
    "#     year, \"_\", szn, \"_mask_6min_pmn.nc\")\n",
    "\n",
    "# # Create dataset object\n",
    "# pmn_ds = Dataset(pmn_filename, \"c\")\n",
    "\n",
    "# # Define lat and lon dims\n",
    "# lon_dim_pm = defDim(pmn_ds, \"lon\", length(mask_ds[\"lon\"][1:subsamp_interval:end]))\n",
    "# lat_dim_pm = defDim(pmn_ds, \"lat\", length(mask_ds[\"lat\"][1:subsamp_interval:end]))\n",
    "\n",
    "# # How to add data to the dims?\n",
    "# lon_dim_pm[:] = mask_ds[\"lon\"][1:subsamp_interval:end]\n",
    "# lat_dim_pm[:] = mask_ds[\"lat\"][1:subsamp_interval:end]\n",
    "\n",
    "# # Define variables pm and pn\n",
    "# pm_var = defVar(pmn_ds, \"pm\", Float32, (\"lon\", \"lat\"))\n",
    "# pn_var = defVar(pmn_ds, \"pn\", Float32, (\"lon\", \"lat\"))\n",
    "\n",
    "# pm_var[:,:] = pm\n",
    "# pn_var[:,:] = pn\n",
    "\n",
    "# println(pmn_ds)\n",
    "\n",
    "# # Close netCDF file\n",
    "# close(pmn_ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export vout as a netCDF file in the same dims as the mask\n",
    "vout_filename = string(output_dir, var_name, \"_\", standard_depth, \"m_\",\n",
    "    year, \"_\", szn, \"_analysis2d.nc\")\n",
    "\n",
    "vout_ds = Dataset(vout_filename, \"c\")\n",
    "\n",
    "# Define lat and lon dims\n",
    "lon_dim_vout = defDim(vout_ds, \"lon\", length(mask_ds[\"lon\"][1:subsamp_interval:end]))\n",
    "lat_dim_vout = defDim(vout_ds, \"lat\", length(mask_ds[\"lat\"][1:subsamp_interval:end]))\n",
    "\n",
    "# How to add data to the dims?\n",
    "lon_dim_vout[:] = mask_ds[\"lon\"][1:subsamp_interval:end]\n",
    "lat_dim_vout[:] = mask_ds[\"lat\"][1:subsamp_interval:end]\n",
    "\n",
    "vout_var = defVar(vout_ds, \"vout\", Float32, (\"lon\", \"lat\"))\n",
    "vout_var[:,:] = vout\n",
    "\n",
    "println(vout_ds)\n",
    "\n",
    "close(vout_ds)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.3",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}