8 Basic xarray - Saving data to netcdf.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Saving data to netCDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "xarray supports saving data to a number of formats, but this tutorial will only cover saving to netCDF files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most common usage pattern for saving data is one where data is read from a netCDF file, manipulated or analysed in some way and then written out again. This is what is assumed in this tutorial.\n",
    "\n",
    "For more information refer to the excellent xarray documentation on [reading and writing files](http://xarray.pydata.org/en/stable/io.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the surface air temperature dataset that has been used previously"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = 'http://dapds00.nci.org.au/thredds/dodsC/rr3/CMIP5/output1/CSIRO-BOM/ACCESS1-3/historical/mon/atmos/Amon/r1i1p1/latest/tas/tas_Amon_ACCESS1-3_historical_r1i1p1_185001-200512.nc'\n",
    "ds = xarray.open_dataset(url)\n",
    "tas = ds.tas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:    (bnds: 2, lat: 145, lon: 192, time: 1872)\n",
       "Coordinates:\n",
       "  * time       (time) datetime64[ns] 1850-01-16T12:00:00 ... 2005-12-16T12:00:00\n",
       "  * lat        (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon        (lon) float64 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1\n",
       "    height     float64 ...\n",
       "Dimensions without coordinates: bnds\n",
       "Data variables:\n",
       "    time_bnds  (time, bnds) datetime64[ns] ...\n",
       "    lat_bnds   (lat, bnds) float64 ...\n",
       "    lon_bnds   (lon, bnds) float64 ...\n",
       "    tas        (time, lat, lon) float32 ...\n",
       "Attributes:\n",
       "    institution:                     CSIRO (Commonwealth Scientific and Indus...\n",
       "    institute_id:                    CSIRO-BOM\n",
       "    experiment_id:                   historical\n",
       "    source:                          ACCESS1-3 2011. Atmosphere: AGCM v1.0 (N...\n",
       "    model_id:                        ACCESS1.3\n",
       "    forcing:                         GHG, Oz, SA, Sl, Vl, BC, OC, (GHG = CO2,...\n",
       "    parent_experiment_id:            piControl\n",
       "    parent_experiment_rip:           r1i1p1\n",
       "    branch_time:                     90945.0\n",
       "    contact:                         The ACCESS wiki: http://wiki.csiro.au/co...\n",
       "    history:                         Fri Apr 13 09:38:12 2012: ncatted -a for...\n",
       "    references:                      See http://wiki.csiro.au/confluence/disp...\n",
       "    initialization_method:           1\n",
       "    physics_version:                 1\n",
       "    tracking_id:                     7f51888d-7daa-45b3-b568-9ce3288b333d\n",
       "    version_number:                  v20120413\n",
       "    product:                         output\n",
       "    experiment:                      historical\n",
       "    frequency:                       mon\n",
       "    creation_date:                   2012-02-05T23:50:03Z\n",
       "    Conventions:                     CF-1.4\n",
       "    project_id:                      CMIP5\n",
       "    table_id:                        Table Amon (27 April 2011) 9c851218e3842...\n",
       "    title:                           ACCESS1-3 model output prepared for CMIP...\n",
       "    parent_experiment:               pre-industrial control\n",
       "    modeling_realm:                  atmos\n",
       "    realization:                     1\n",
       "    cmor_version:                    2.8.0\n",
       "    DODS_EXTRA.Unlimited_Dimension:  time"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because xarray saves all the metadata and coordinates in a dataset object it is very simple to write it back out to a netCDF file. To save a dataset to a local file call the `to_netcdf` method and supply a filename"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.to_netcdf('ds.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that is it. The saved file can be read back into a new variable, `ds_local` and inspection of the metadata suggests it is the same. It is impossible to say for certain just by inspecting the metadata, luckily ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_local = xarray.open_dataset('ds.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:    (bnds: 2, lat: 145, lon: 192, time: 1872)\n",
       "Coordinates:\n",
       "  * time       (time) datetime64[ns] 1850-01-16T12:00:00 ... 2005-12-16T12:00:00\n",
       "  * lat        (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon        (lon) float64 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1\n",
       "    height     float64 ...\n",
       "Dimensions without coordinates: bnds\n",
       "Data variables:\n",
       "    time_bnds  (time, bnds) datetime64[ns] ...\n",
       "    lat_bnds   (lat, bnds) float64 ...\n",
       "    lon_bnds   (lon, bnds) float64 ...\n",
       "    tas        (time, lat, lon) float32 ...\n",
       "Attributes:\n",
       "    institution:                     CSIRO (Commonwealth Scientific and Indus...\n",
       "    institute_id:                    CSIRO-BOM\n",
       "    experiment_id:                   historical\n",
       "    source:                          ACCESS1-3 2011. Atmosphere: AGCM v1.0 (N...\n",
       "    model_id:                        ACCESS1.3\n",
       "    forcing:                         GHG, Oz, SA, Sl, Vl, BC, OC, (GHG = CO2,...\n",
       "    parent_experiment_id:            piControl\n",
       "    parent_experiment_rip:           r1i1p1\n",
       "    branch_time:                     90945.0\n",
       "    contact:                         The ACCESS wiki: http://wiki.csiro.au/co...\n",
       "    history:                         Fri Apr 13 09:38:12 2012: ncatted -a for...\n",
       "    references:                      See http://wiki.csiro.au/confluence/disp...\n",
       "    initialization_method:           1\n",
       "    physics_version:                 1\n",
       "    tracking_id:                     7f51888d-7daa-45b3-b568-9ce3288b333d\n",
       "    version_number:                  v20120413\n",
       "    product:                         output\n",
       "    experiment:                      historical\n",
       "    frequency:                       mon\n",
       "    creation_date:                   2012-02-05T23:50:03Z\n",
       "    Conventions:                     CF-1.4\n",
       "    project_id:                      CMIP5\n",
       "    table_id:                        Table Amon (27 April 2011) 9c851218e3842...\n",
       "    title:                           ACCESS1-3 model output prepared for CMIP...\n",
       "    parent_experiment:               pre-industrial control\n",
       "    modeling_realm:                  atmos\n",
       "    realization:                     1\n",
       "    cmor_version:                    2.8.0\n",
       "    DODS_EXTRA.Unlimited_Dimension:  time"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds_local"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is an operator, `identical`, that will compare datasets and return `True` if they are identical. Running that shows that yes, xarray has saved the dataset read from an opendap URL to a local file and `ds` and `ds_local` are identical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds.identical(ds_local)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same `to_netcdf` operator will work for a `DataArray`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "tas.to_netcdf('tas.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When this is read back in and printed out none of the metadata from the original dataset, `ds`, is retained as it did not exist in the `tas` `DataArray`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (lat: 145, lon: 192, time: 1872)\n",
       "Coordinates:\n",
       "  * time     (time) datetime64[ns] 1850-01-16T12:00:00 ... 2005-12-16T12:00:00\n",
       "  * lat      (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon      (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n",
       "    height   float64 ...\n",
       "Data variables:\n",
       "    tas      (time, lat, lon) float32 ..."
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tas_local = xarray.open_dataset('tas.nc')\n",
    "tas_local"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the `identical` method shows the newly created `tas` `Dataset` is not identical to the `DataArray` `tas`. If the comparison is between `tas` and the `DataArray` `tas` inside the `tas_local` `Dataset` it shows they are identical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tas.identical(tas_local)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tas.identical(tas_local.tas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `DataArray` cannot be written directly to a netCDF file. When `to_netcdf` is called on a `DataArray` there is an implicit `to_dataset` operation then `to_netcdf` is called on the resulting `DataSet`. So the result of `tas.to_dataset()` is identical to `tas_local`, as this is what was done implicitly before `to_netcdf`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tas.to_dataset().identical(tas_local)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To save new datasets, say the results of calculations, is the same process. As an example from the surface air temp data, `tas`, a 30 year climatology is calculated and saved in the `tas_clim` variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "tas_clim = tas.sel(time=slice('1960-01','1989-12')).mean(dim='time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'tas' (lat: 145, lon: 192)>\n",
       "array([[223.5951 , 223.5951 , 223.5951 , ..., 223.58675, 223.58675, 223.58675],\n",
       "       [227.05313, 227.00461, 226.95636, ..., 227.20116, 227.1515 , 227.10295],\n",
       "       [228.13997, 228.0243 , 227.9114 , ..., 228.49681, 228.37773, 228.25885],\n",
       "       ...,\n",
       "       [255.42877, 255.46603, 255.50294, ..., 255.29857, 255.34636, 255.38702],\n",
       "       [254.99307, 255.0079 , 255.0223 , ..., 254.94296, 254.95801, 254.97554],\n",
       "       [254.3218 , 254.3218 , 254.3218 , ..., 254.3218 , 254.3218 , 254.3218 ]],\n",
       "      dtype=float32)\n",
       "Coordinates:\n",
       "  * lat      (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon      (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n",
       "    height   float64 1.5"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tas_clim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now this is saved to the netCDF file `tas_clim.nc`. Note that this is calling `to_dataset` explicitly so that a new name can be specified for the `DataArray` when it is converted into a `Dataset`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "tas_clim.to_dataset(name='tas_climatology').to_netcdf('tas_clim.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open the data again, overwriting the `tas_clim` variable in the process, and the data variable is now named `tas_climatology` as expected. If the name had not been defined it would have used the name from the original `DataArray`, namely `tas`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "tas_clim = xarray.open_dataset('tas_clim.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:          (lat: 145, lon: 192)\n",
       "Coordinates:\n",
       "  * lat              (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 87.5 88.75 90.0\n",
       "  * lon              (lon) float64 0.0 1.875 3.75 5.625 ... 354.4 356.2 358.1\n",
       "    height           float64 ...\n",
       "Data variables:\n",
       "    tas_climatology  (lat, lon) float32 ..."
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tas_clim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When xarray opens a netCDF file some of the metadata is saved in a special `dict` called `encoding`. This is used when xarray writes a `Dataset` back out to a netCDF file, it encodes the data before doing so. There is an `encoding` attribute for a `Dataset` and all variables and coordinates have their own `encoding` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'unlimited_dims': {'time'},\n",
       " 'source': 'http://dapds00.nci.org.au/thredds/dodsC/rr3/CMIP5/output1/CSIRO-BOM/ACCESS1-3/historical/mon/atmos/Amon/r1i1p1/latest/tas/tas_Amon_ACCESS1-3_historical_r1i1p1_185001-200512.nc'}"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds.encoding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'source': 'http://dapds00.nci.org.au/thredds/dodsC/rr3/CMIP5/output1/CSIRO-BOM/ACCESS1-3/historical/mon/atmos/Amon/r1i1p1/latest/tas/tas_Amon_ACCESS1-3_historical_r1i1p1_185001-200512.nc',\n",
       " 'original_shape': (1872, 145, 192),\n",
       " 'dtype': dtype('float32'),\n",
       " 'missing_value': 1e+20,\n",
       " '_FillValue': 1e+20,\n",
       " 'coordinates': 'height'}"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tas.encoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `encoding` `dict` for the `time` coordinate is particularly interesting, as this is where the `calendar` and `units` attributes from the netCDF file are stored. When the netCDF file is opened these attributes are used to convert the time variable, which is just an array of numbers, into python `DateTime` objects. When the data is written back to a file this step must be done in reverse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'source': 'http://dapds00.nci.org.au/thredds/dodsC/rr3/CMIP5/output1/CSIRO-BOM/ACCESS1-3/historical/mon/atmos/Amon/r1i1p1/latest/tas/tas_Amon_ACCESS1-3_historical_r1i1p1_185001-200512.nc',\n",
       " 'original_shape': (1872,),\n",
       " 'dtype': dtype('float64'),\n",
       " 'units': 'days since 0001-01-01',\n",
       " 'calendar': 'proleptic_gregorian'}"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds.time.encoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "netCDF4 files support lossless in-file compression. This can reduce the size of a netCDF file considerably, but it can still be used in the same way as an uncompressed netCDF file. See the CMS Wiki for [more details in netCDF compression](http://climate-cms.wikis.unsw.edu.au/NetCDF_Compression_Tools).\n",
    "\n",
    "Enabling in-file compression when saving data to a netCDF file is done by adding the appropriate key/value pairs to the `encoding`. The `encoding` can be updated in the object before writing to a file, or an `encoding` argument added to the `to_netcdf` method. The latter is shown below for the the `tas` `DataArray`. It is sufficient to just specify `zlib`, which would use the default compression level, `complevel`, of 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "tas.to_netcdf('tas_compressed.nc', encoding=tas.encoding.update({'zlib': True, 'complevel': 4}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A directory listing shows that the compressed file is 56% the size of the uncompressed data, a 1.8 compression factor. 3D ocean data will typically have compression ratios of 3-4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r-- 1 aph502 v45 111M Jul 25 16:02 tas_compressed.nc\r\n",
      "-rw-r--r-- 1 aph502 v45 199M Jul 25 15:56 tas.nc\r\n"
     ]
    }
   ],
   "source": [
    "!ls -lh tas.nc tas_compressed.nc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reading the compressed data back in shows the compression was lossless and the data is unchanged"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tas.identical(xarray.open_dataset('tas_compressed.nc').tas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same `encoding` option works for `Dataset` objects also"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.to_netcdf('ds_compressed.nc', encoding=ds.encoding.update({'zlib': True, 'complevel': 4}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r-- 1 aph502 v45 113M Jul 25 16:03 ds_compressed.nc\r\n",
      "-rw-r--r-- 1 aph502 v45 200M Jul 25 15:55 ds.nc\r\n"
     ]
    }
   ],
   "source": [
    "!ls -lh ds.nc ds_compressed.nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:analysis3-19.07]",
   "language": "python",
   "name": "conda-env-analysis3-19.07-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}