03_Postprocessing_Combine_Netcdfs_per_season.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Script to combine all the separate NetCDF files for one species per season\n",
    "\n",
    "* The script loops over the variable names in _varname_\n",
    "* creates a new directory /_varname_/netcdf_season/\n",
    "* and combines the data of the netcdfs in /_varname_/netcdf_all into 4 netcdfs (one per season): varname_StartyearEndyear_season.nc\n",
    "\n",
    "with as season a number between 1 and 4:\n",
    "\n",
    "* 1: months 1-3\n",
    "* 2: months 4-6\n",
    "* 3: months 7-9\n",
    "* 4: months 10-12\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Julia Version 0.6.4\n",
      "Commit 9d11f62bcb (2018-07-09 19:09 UTC)\n",
      "Platform Info:\n",
      "  OS: Linux (x86_64-pc-linux-gnu)\n",
      "  CPU: Intel(R) Xeon(R) CPU E5-2630L 0 @ 2.00GHz\n",
      "  WORD_SIZE: 64\n",
      "  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge MAX_THREADS=16)\n",
      "  LAPACK: libopenblas64_\n",
      "  LIBM: libopenlibm\n",
      "  LLVM: libLLVM-3.9.1 (ORCJIT, sandybridge)\n"
     ]
    }
   ],
   "source": [
    "versioninfo()\n",
    "print(join([\"$p: $v\\n\" for (p,v) in Pkg.installed()]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/data/20180306_OOPS\""
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Libz: 0.4.0\n",
      "AxisAlgorithms: 0.3.0\n",
      "OffsetArrays: 0.6.0\n",
      "HTTP: 0.6.12\n",
      "CondaBinDeps: 0.1.0\n",
      "OceanPlot: 0.0.0-\n",
      "Nullables: 0.0.6\n",
      "ZMQ: 0.6.3\n",
      "DataStructures: 0.8.4\n",
      "Compat: 1.0.0\n",
      "EzXML: 0.7.1\n",
      "ShowItLikeYouBuildIt: 0.2.0\n",
      "SpecialFunctions: 0.6.0\n",
      "Blosc: 0.5.0\n",
      "ZipFile: 0.6.0\n",
      "FixedPointNumbers: 0.4.6\n",
      "SHA: 0.5.7\n",
      "Missings: 0.2.10\n",
      "MAT: 0.4.0\n",
      "HDF5: 0.9.3\n",
      "ColorTypes: 0.6.7\n",
      "BufferedStreams: 0.4.1\n",
      "MbedTLS: 0.5.12\n",
      "SortingAlgorithms: 0.2.1\n",
      "Conda: 1.0.0\n",
      "PyCall: 1.17.1\n",
      "WoodburyMatrices: 0.3.0\n",
      "JSON: 0.17.2\n",
      "StatsBase: 0.23.1\n",
      "IJulia: 1.9.1\n",
      "PyPlot: 2.6.0\n",
      "BinDeps: 0.8.8\n",
      "Parameters: 0.9.1\n",
      "DIVAnd: 2.0.0+\n",
      "Mustache: 0.3.3\n",
      "Reexport: 0.1.0\n",
      "CMakeWrapper: 0.1.0\n",
      "URIParser: 0.3.1\n",
      "Interpolations: 0.8.0\n",
      "BinaryProvider: 0.3.3\n",
      "IniFile: 0.4.0\n",
      "LaTeXStrings: 1.0.0\n",
      "DataArrays: 0.7.0\n",
      "VersionParsing: 1.1.1\n",
      "Colors: 0.8.2\n",
      "MacroTools: 0.4.2\n",
      "NCDatasets: 0.3.0\n",
      "Ratios: 0.2.0\n"
     ]
    }
   ],
   "source": [
    "pwd()\n",
    "cd(\"//data/20180306_OOPS/\")\n",
    "pwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "using NCDatasets\n",
    "using DIVAnd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1 create netcdf from template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Adjusted from the output of the function 'ncgen'\n",
    "# ncgen(\"test.nc\", newfname = \"test3.nc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "create_netcdf_season_template (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function create_netcdf_season_template(varname,yearrange,season, varncsDir)\n",
    "    # Adjusted from the output of the function 'ncgen'\n",
    "    # ncgen(\"test.nc\", newfname = \"test3.nc\")\n",
    "    \n",
    "    ds = Dataset(joinpath(varncsDir,\"$(replace(varname,\"log_\",\"\"))_$(yearrange[1])$(yearrange[end])_$(season).nc\"),\"c\")\n",
    "    println(\"creating netcdf: $(replace(varname,\"log_\",\"\"))_$(yearrange[1])$(yearrange[end])_$(season).nc\")   \n",
    "    \n",
    "    # Dimensions\n",
    "\n",
    "    ds.dim[\"lon\"] = 951\n",
    "    ds.dim[\"lat\"] = 401\n",
    "    ds.dim[\"time\"] = length(yearrange)\n",
    "    ds.dim[\"nv\"] = 2\n",
    "    #ds.dim[\"season\"] = 4\n",
    "    # ds.dim[\"year\"] = length(yearrange)\n",
    "\n",
    "    # Declare variables\n",
    "\n",
    "    nclon = defVar(ds,\"lon\", Float64, (\"lon\",)) \n",
    "    nclon.attrib[\"units\"] = \"degrees_east\"\n",
    "    nclon.attrib[\"standard_name\"] = \"longitude\"\n",
    "    nclon.attrib[\"long_name\"] = \"longitude\"\n",
    "\n",
    "    nclat = defVar(ds,\"lat\", Float64, (\"lat\",)) \n",
    "    nclat.attrib[\"units\"] = \"degrees_north\"\n",
    "    nclat.attrib[\"standard_name\"] = \"latitude\"\n",
    "    nclat.attrib[\"long_name\"] = \"latitude\"\n",
    "\n",
    "    nctime = defVar(ds,\"time\", Float64, (\"time\",)) \n",
    "    nctime.attrib[\"units\"] = \"days since 1900-01-01 00:00:00\"\n",
    "    nctime.attrib[\"standard_name\"] = \"time\"\n",
    "    nctime.attrib[\"long_name\"] = \"time\"\n",
    "    nctime.attrib[\"calendar\"] = \"standard\"\n",
    "    nctime.attrib[\"climatology\"] = \"climatology_bounds\"\n",
    "\n",
    "    ncclimatology_bounds = defVar(ds,\"climatology_bounds\", Float64, (\"nv\", \"time\",)) \n",
    "    ncclimatology_bounds.attrib[\"units\"] = \"days since 1900-01-01 00:00:00\"\n",
    "\n",
    "    # ncseason = defVar(ds,\"season\", Int, (\"season\",))\n",
    "    # ncseason.attrib[\"units\"] = \"season\"\n",
    "    # ncseason.attrib[\"standard_name\"] = \"season\"\n",
    "    # ncseason.attrib[\"long_name\"] = \"season number (1: months 1-3; 2: months 4-6; 3: months 7-9; 4: months 10-12)\"\n",
    "    \n",
    "    # ncyear = defVar(ds,\"year\", Int, (\"year\",))\n",
    "    # ncyear.attrib[\"units\"] = \"year\"\n",
    "    # ncyear.attrib[\"standard_name\"] = \"year\"\n",
    "    # ncyear.attrib[\"long_name\"] = \"year\"\n",
    "    # ncyear.attrib[\"calendar\"] = \"standard\"\n",
    "    \n",
    "    nclog = defVar(ds,\"$(varname)\", Float32, (\"lon\", \"lat\", \"time\")) # removed \"time\", \"season\"\n",
    "    nclog.attrib[\"units\"] = \"logab\"\n",
    "    nclog.attrib[\"standard_name\"] = \"$(varname)\"\n",
    "    nclog.attrib[\"long_name\"] = \"$(replace(varname,'_',' '))\"\n",
    "    nclog.attrib[\"_FillValue\"] = Float32(9.96921e36)\n",
    "    nclog.attrib[\"missing_value\"] = Float32(9.96921e36)\n",
    "    nclog.attrib[\"cell_methods\"] = \"time: mean within years time: mean over years\"\n",
    "\n",
    "    nclog_L1 = defVar(ds,\"$(varname)_L1\", Float32, (\"lon\", \"lat\", \"time\")) # removed \"time\", \"season\"\n",
    "    nclog_L1.attrib[\"units\"] = \"logab\"\n",
    "    nclog_L1.attrib[\"standard_name\"] = \"$(varname)_L1\"\n",
    "    nclog_L1.attrib[\"long_name\"] = \"$(replace(varname,'_',' ')) masked using relative error threshold 0.3\"\n",
    "    nclog_L1.attrib[\"_FillValue\"] = Float32(9.96921e36)\n",
    "    nclog_L1.attrib[\"missing_value\"] = Float32(9.96921e36)\n",
    "\n",
    "    nclog_L2 = defVar(ds,\"$(varname)_L2\", Float32, (\"lon\", \"lat\", \"time\")) # removed \"time\", \"season\"\n",
    "    nclog_L2.attrib[\"units\"] = \"logab\"\n",
    "    nclog_L2.attrib[\"standard_name\"] = \"$(varname)_L2\"\n",
    "    nclog_L2.attrib[\"long_name\"] = \"$(replace(varname,'_',' ')) masked using relative error threshold 0.5\"\n",
    "    nclog_L2.attrib[\"_FillValue\"] = Float32(9.96921e36)\n",
    "    nclog_L2.attrib[\"missing_value\"] = Float32(9.96921e36)\n",
    "\n",
    "    nclog_relerr = defVar(ds,\"$(varname)_relerr\", Float32, (\"lon\", \"lat\", \"time\")) # removed \"time\", \"season\"\n",
    "    nclog_relerr.attrib[\"units\"] = \"1\"\n",
    "    nclog_relerr.attrib[\"long_name\"] = \"Relative error of $(replace(varname,'_',' '))\"\n",
    "    nclog_relerr.attrib[\"valid_min\"] = Float32(0.0)\n",
    "    nclog_relerr.attrib[\"valid_max\"] = Float32(1.0)\n",
    "    nclog_relerr.attrib[\"_FillValue\"] = Float32(9.96921e36)\n",
    "    nclog_relerr.attrib[\"missing_value\"] = Float32(9.96921e36)\n",
    "\n",
    "    # Global attributes\n",
    "\n",
    "    ds.attrib[\"Conventions\"] = \"CF-1.6\"\n",
    "    ds.attrib[\"title\"] = \"DIVA 3D analysis of $(replace(varname,'_',' '))\"\n",
    "    ds.attrib[\"file_name\"] = \"$(varname).nc\"\n",
    "    ds.attrib[\"product_id\"] = \"16342e64-c6e3-11e8-320c-055e34a07963\"\n",
    "    ds.attrib[\"date\"] = \"2018-10-03T10:05:26\"\n",
    "    ds.attrib[\"project\"] = \"EMODnet Biology\"\n",
    "    ds.attrib[\"institution_urn\"] = \"SDN::EDMO::422\"\n",
    "    ds.attrib[\"production\"] = \"Flanders Marine Institute\"\n",
    "    ds.attrib[\"Author_e-mail\"] = \"Lennert Schepers <lennert.schepers@vliz.be>\"\n",
    "    ds.attrib[\"source\"] = \"SAHFOS / Marine Biological Association (UK)\"\n",
    "    ds.attrib[\"comment\"] = \"...\"\n",
    "    ds.attrib[\"product_version\"] = \"2.0\"\n",
    "\n",
    "    close(ds)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### populate netcdf  - per season"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@show t1 = now()\n",
    "\n",
    "yearrange = 1958:2016\n",
    "for varname = [\"log_tem_lon\",\"log_cal_hel\", \"log_tot_lar\",\"log_tot_sma\",\n",
    "               \"log_oit_spp\", \"log_cal_fin\",\"ratio_large_to_small\",\n",
    "               \"log_aca_spp\", \"log_met_luc\", \"log_chli\"]\n",
    "    \n",
    "    varDir = joinpath(pwd(),varname)             # species directory\n",
    "    varncDir = joinpath(varDir,\"netcdf_all\")     #  subdirectory with all separate netcdfs\n",
    "    varncsDir = joinpath(varDir,\"netcdf_season\") #  subdirectory for combined netcdfs\n",
    "    mkpath(varncsDir)\n",
    "    \n",
    "    for season = 1:4\n",
    "                # read data\n",
    "                print(\"season: $(season)  \")\n",
    "                # define monthstring\n",
    "                if season == 1 \n",
    "                    season_str = \"0103\"\n",
    "                elseif season == 2\n",
    "                    season_str = \"0406\"\n",
    "                elseif season == 3\n",
    "                    season_str = \"0709\"\n",
    "                elseif season == 4\n",
    "                    season_str = \"1012\"\n",
    "                end\n",
    "    \n",
    "    \n",
    "        i = 1\n",
    "        # create template\n",
    "        println(\"creating netcdf template ...\")\n",
    "        create_netcdf_season_template(varname, yearrange, season, varncsDir)\n",
    "        \n",
    "\n",
    "        # append data to the NetCDF\n",
    "        Dataset(joinpath(varncsDir,\"$(replace(varname,\"log_\",\"\"))_$(yearrange[1])$(yearrange[end])_$(season).nc\"),\"a\") do ds \n",
    "\n",
    "            for year = yearrange\n",
    "                print(\"\\nyear: $(year) \")\n",
    "                \n",
    "                fname = joinpath(varncDir,\"$(replace(varname,\"log_\",\"\"))_$(year)$(year)_$(season_str).nc\")\n",
    "                \n",
    "                # variables only needed 1st time\n",
    "                if i==1\n",
    "                    Dataset(fname,\"r\") do nc\n",
    "                        # take values of netcdf\n",
    "                        ds[\"lat\"][:] = nc[\"lat\"][:]\n",
    "                        ds[\"lon\"][:] = nc[\"lon\"][:]\n",
    "                        \n",
    "                        # ds[\"season\"][:] = 1:4\n",
    "                        # ds[\"year\"][:] = yearrange\n",
    "                        \n",
    "                        ds.attrib[\"product_id\"] = nc.attrib[\"product_id\"]\n",
    "                        \n",
    "                    end\n",
    "                end # end if\n",
    "                \n",
    "                # for each netcdf\n",
    "                Dataset(fname,\"r\") do nc\n",
    "\n",
    "                    # take values of slice\n",
    "                    ds[\"time\"][i] = nc[\"time\"][1] ## HOW IS TIME CONVERTED TO FLOAT?\n",
    "                    ds[\"climatology_bounds\"][:,i] = nc[\"climatology_bounds\"].var[:,1]\n",
    "                    \n",
    "                    # \"lon\", \"lat\", \"time\", \"season\", \"year\"\n",
    "                    ds[\"$(varname)\"][:,:, i] = nc[\"$(varname)\"][:,:,1,1]\n",
    "                    ds[\"$(varname)_L1\"][:,:, i] = nc[\"$(varname)_L1\"][:,:,1,1]\n",
    "                    ds[\"$(varname)_L2\"][:,:, i] = nc[\"$(varname)_L2\"][:,:,1,1]\n",
    "                    ds[\"$(varname)_relerr\"][:,:, i] = nc[\"$(varname)_relerr\"][:,:,1,1]\n",
    "                    \n",
    "                end # end single netcdf\n",
    "                \n",
    "\n",
    "            i = i + 1  \n",
    "         \n",
    "            end # year        \n",
    "        end # ds closed\n",
    "   end # season\n",
    "  \n",
    "end # varname\n",
    "@show t2=now()\n",
    "@show Dates.canonicalize(Dates.CompoundPeriod(t2-t1))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.4",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}