data_files.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "# Choosing the right data file format for numerical computing\n",
    "\n",
    "This notebook will go over the pros and cons of various data file formats common in numerical python workflows. It'll cover various concerns when storing data on-disk and how popular formats address these challenges; the what and why of these formats."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "## Who am I?\n",
    "\n",
    "* David Hoese <sub>pronounced like Haze</sub>\n",
    "* Software Developer (Sr. Instrumentation Technologist)<br>\n",
    "  Space Science and Engineering Center (SSEC)<br>\n",
    "  University of Wisconsin - Madison\n",
    "* Satpy and Vispy developer\n",
    "* @djhoese on Twitter and GitHub\n",
    "\n",
    "This notebook:\n",
    "https://github.com/djhoese/data-file-formats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "### Why write data to disk?\n",
    "\n",
    "1. State or Cache\n",
    "   * a long running application needs to start where it left off\n",
    "   * reuse calculation in future executions\n",
    "   * user settings/preferences are saved for later use\n",
    "2. Data Archival/Distribution/Sharing\n",
    "   * results are shared with other people\n",
    "   * results are shared with other software"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "### What kind of files are we working with?\n",
    "\n",
    "* Plain text or binary file formats\n",
    "* Primarily numeric data\n",
    "* Optionally add custom metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "### What are we **NOT** talking about?\n",
    "\n",
    "* Databases\n",
    "* Python's pickle format\n",
    "* Storing user application settings\n",
    "* **Custom** binary formats (no community/organization support)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why does file format matter?\n",
    "\n",
    "There are many different formats that we can use in Python and we have\n",
    "different ways to access them. Each one comes with various advantages and\n",
    "disadvantages. When it comes to choosing a file format for a task we should\n",
    "be concerned with a few key things:\n",
    "\n",
    "* **Write speeds**<br>\n",
    "  How long does it take to get data from it's in-memory format to disk?\n",
    "* **Read speeds**<br>\n",
    "  How long does it take to get data from disk to a usable in-memory format?\n",
    "* **On-disk size**<br>\n",
    "  How big are the files?\n",
    "* **In-memory size**<br>\n",
    "  How much memory does it take to open the file and get data out? Are there multiple copies of the data in memory?\n",
    "* **Flexibility/Usefulness**<br>\n",
    "  Can I do anything else with the format? Store metadata? Can I use it for this other use case? Can I lazily load data? Can I access data remotely? Can I compress the data? How much data can I store in one file?\n",
    "* **Format Availability**<br>\n",
    "  Is the format only readable by Python software? Do my users need to learn something new to read it? If the format isn't \"built in\" to a programming environment, can I easily install the necessary libraries?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Warmup - Plain Text\n",
    "\n",
    "Let's start with an example that uses plain text to write the last time that\n",
    "the function was called. We first need to import a few modules and then\n",
    "create our function to append to our data file called `last_time.txt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "DATA_FILE = 'last_time.txt'\n",
    "\n",
    "def write_time():\n",
    "    with open(DATA_FILE, 'a') as time_file:\n",
    "        now = datetime.utcnow()\n",
    "        data_line = \"{:f}\\n\".format(now.timestamp())\n",
    "        time_file.write(data_line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our function we get the current time in\n",
    "[UTC](https://en.wikipedia.org/wiki/Coordinated_Universal_Time). We convert\n",
    "it to epoch time (a.k.a. POSIX time), seconds since `1970-01-01 00:00:00`, and\n",
    "write the string representation to the file as a single line.\n",
    "\n",
    "Let's add a couple times to the file by running the `write_time` a couple\n",
    "times below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "write_time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can use the command `head` to print out the first couple lines of the file's contents:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head $DATA_FILE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've written data to the file, now let's read it back out. The below function\n",
    "will read the content's of the file in to a list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_times():\n",
    "    exec_times = []\n",
    "    with open(DATA_FILE, 'r') as time_file:\n",
    "        for line in time_file:\n",
    "            time_val = float(line.strip())\n",
    "            exec_times.append(time_val)\n",
    "    return exec_times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "read_times()[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code above gives us a simple example of saving data from software to a\n",
    "file on disk. We wrote a single value at a time and accumulate more\n",
    "information as time went on. We were able to read these data back\n",
    "in to python at a later time. Could we have done anything differently?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plain Text - Pros/Cons\n",
    "\n",
    "Let's take the above concerns and look at our text file from before and the\n",
    "code we used to access it.\n",
    "\n",
    "<details>\n",
    "<summary>Pros</summary>\n",
    "    \n",
    "* Human readable\n",
    "* Simple code (no external libraries)\n",
    "* Easily usable by other languages/tools\n",
    "* Could read one value at a time (but we didn't)\n",
    "    \n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Cons</summary>\n",
    "\n",
    "* Have to convert data to/from string and float (slow)\n",
    "* Representing each 8 byte float (64-bit) as ~17 ASCII bytes\n",
    "* Unknown precision of data values, how many decimal points?\n",
    "* Don't know how many elements until it is fully read\n",
    "* Can't easily seek to a specific index/element\n",
    "* Code: Read as a list instead of a numpy array and used a python for loop (potentially slow)\n",
    "    \n",
    "</details>\n",
    "<br>\n",
    "\n",
    "And here's what a single value of our text file looks like on disk (8-bit ASCII character):\n",
    "\n",
    "<table>\n",
    "    <tr>\n",
    "        <th>Byte Offset</th>\n",
    "        <th>Value</th>\n",
    "    </tr>\n",
    "    <tr><td>0</td><td>1</td></tr>\n",
    "    <tr><td>1</td><td>5</td></tr>\n",
    "    <tr><td>2</td><td>6</td></tr>\n",
    "    <tr><td>3</td><td>8</td></tr>\n",
    "    <tr><td>4</td><td>6</td></tr>\n",
    "    <tr><td>5</td><td>6</td></tr>\n",
    "    <tr><td>6</td><td>7</td></tr>\n",
    "    <tr><td>7</td><td>8</td></tr>\n",
    "    <tr><td>8</td><td>3</td></tr>\n",
    "    <tr><td>9</td><td>2</td></tr>\n",
    "    <tr><td>10</td><td>.</td></tr>\n",
    "    <tr><td>11</td><td>8</td></tr>\n",
    "    <tr><td>12</td><td>2</td></tr>\n",
    "    <tr><td>13</td><td>7</td></tr>\n",
    "    <tr><td>14</td><td>4</td></tr>\n",
    "    <tr><td>15</td><td>6</td></tr>\n",
    "    <tr><td>16</td><td>8</td></tr>\n",
    "    <tr><td>17</td><td>\\n</td></tr>\n",
    "    </table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can perform a quick timing to see how long it takes to read the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "txt_time_write = %timeit -o -n 10000 -r 1 write_time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "txt_time_read = %timeit -o -n 100 -r 1 read_times()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NumPy library also provides a function for loading data from text files.\n",
    "Let's try it and see how it compares."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "txt_time_read_np = %timeit -o -n 100 -r 1 np.loadtxt(DATA_FILE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems that in this specific case and with all of the extra checks numpy\n",
    "performs, it actually takes longer to read the data with numpy. When it comes\n",
    "to simple human-readable formats, we couldn't have gone much simpler.\n",
    "\n",
    "The remainder of this document will go through different use cases and\n",
    "the file formats that we have as options. We'll apply this type of\n",
    "performance analysis to our format choices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Flat Binary\n",
    "\n",
    "When human-readability isn't necessary, another option for storing simple\n",
    "data structures is a flat binary file. A flat binary file consists of the\n",
    "raw data values stored contiguously as a flat array. Let's rewrite our code\n",
    "from above to write to a flat binary file using the numpy library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "BIN_DATA_FILE = 'last_time.dat'\n",
    "\n",
    "def write_time_binary():\n",
    "    with open(BIN_DATA_FILE, 'a') as time_file:\n",
    "        now = datetime.utcnow()\n",
    "        np.array([now.timestamp()], dtype=np.float64).tofile(time_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_times_binary():\n",
    "     return np.fromfile(BIN_DATA_FILE, dtype=np.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_time_write = %timeit -o -n 100000 -r 1 write_time_binary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_time_read = %timeit -o -n 100 -r 7 read_times_binary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform Tip: Memory Maps\n",
    "\n",
    "By default, when reading a file from disk we have to transfer data from\n",
    "the hard disk to system memory. That means that when creating something\n",
    "like a numpy array from a binary data file we are transferring **all** of the\n",
    "file's contents from disk in to memory when we might not use it right away; a very slow operation. There is a\n",
    "lazier, generally more efficient, method called memory mapping (a.k.a. mmap in\n",
    "C, memmap by numpy). By creating a memory map, we allocate the virtual memory\n",
    "space for our data, but the operating system won't load the data from disk\n",
    "until we need it. Memory mapping avoids extra copies of file data in memory, works very well with random accesses to data in a file, can be cached more efficiently, shared between processes more effectively, and is generally your best option for reading large files that are difficult to hold in memory at a single time.\n",
    "\n",
    "Further reading:\n",
    "\n",
    "* [Stackoverflow answer discussing memory maps](https://stackoverflow.com/a/6383253/433202)\n",
    "* [Memory-mapped File - Wikipedia](https://en.wikipedia.org/wiki/Memory-mapped_file)\n",
    "\n",
    "Going back to the above binary file usage..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_times_binary_memmap():\n",
    "    return np.memmap(BIN_DATA_FILE, mode='r', dtype=np.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_time_read_mmap = %timeit -o -n 100 -r 7 read_times_binary_memmap()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_time_read.average / bin_time_read_mmap.average"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keep in mind that memory mapping isn't **loading** the data in to memory so this isn't technically a fair comparison. However, as we use the memory mapped array object we should see better performance than a traditional read of the file.\n",
    "\n",
    "Note that we could also use memory maps to write data. This is most beneficial if we are writing to random locations in the file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "### Flat Binary - Pros/Cons\n",
    "\n",
    "<details>\n",
    "    <summary>Pros</summary>\n",
    "    \n",
    "* Simple code\n",
    "* Readable by any programming language (*see Cons)\n",
    "* Minimum on-disk size without compression\n",
    "* Fast reading and memory mappable\n",
    "* Supports N-dimensional data\n",
    "* Subsettable\n",
    "\n",
    "</details>\n",
    "    \n",
    "<details>\n",
    "    <summary>Cons</summary>\n",
    "    \n",
    "* Not human readable\n",
    "* No shape, data type, or byte order information stored\n",
    "* Platform dependent (not shareable)\n",
    "\n",
    "Note from numpy docs:\n",
    "    \n",
    "> Do not rely on the combination of tofile and fromfile for data storage, as the binary files generated are are not platform independent. In particular, no byte-order or data-type information is saved. Data can be stored in the platform independent .npy format using save and load instead.\n",
    "    \n",
    "</details>\n",
    "\n",
    "Storage layout (64-bit float):\n",
    "\n",
    "<table>\n",
    "    <tr>\n",
    "        <th>Byte Offset</th>\n",
    "        <th>Value</th>\n",
    "    </tr>\n",
    "    <tr><td>0</td><td>1568667832.827468</td></tr>\n",
    "    <tr><td>8</td><td>1568667832.827796</td></tr>\n",
    "    <tr><td>16</td><td>1568667832.827835</td></tr>\n",
    "    </table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bonus - Numpy's .npy format\n",
    "\n",
    "As mention in the cons above, a flat binary file *only* has the data and no\n",
    "other information. This means we have to keep track of this information\n",
    "ourselves. To help with this numpy provides a `.npy` format which stores this\n",
    "information inside the file alongside the binary data. Here's a quick example\n",
    "to create a `.npy` file using\n",
    "[`np.save`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.save.html#numpy.save)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('last_time.npy', np.array([1., 2., 3.]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we are ready to read the data back we can use\n",
    "[`np.load`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.load.html#numpy.load).\n",
    "This method gives us the option to open the data as a memory map, giving us\n",
    "the space and speed advantages of a flat binary while avoiding the format\n",
    "metadata issues. Keep in mind that this format is only readable by numpy and\n",
    "would require an additional library in any other language."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.load('last_time.npy', mmap_mode='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Comma-Separated Values (CSV)\n",
    "\n",
    "So far we've dealt with a single stream (a.k.a. field or variable) of data,\n",
    "but what if we need to store more? When it comes to multiple 1-dimensional\n",
    "variables one of the more common solutions is a Comma-Separate Values (CSV)\n",
    "file. We could use numpy again, but instead we'll use the `pandas` library\n",
    "for its more powerful handling of tabular data like we would store in a CSV.\n",
    "\n",
    "We'll start by loading some example data used by the seaborn python package.\n",
    "Their example data is stored as a\n",
    "[CSV file on GitHub](https://raw.githubusercontent.com/mwaskom/seaborn-data/master/iris.csv).\n",
    "We use the pandas\n",
    "[`read_csv`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.read_csv.html)\n",
    "function. This function has a lot of options, but we'll use the defaults for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "seaborn_iris_url = 'https://raw.githubusercontent.com/mwaskom/seaborn-data/master/iris.csv'\n",
    "data = pd.read_csv(seaborn_iris_url)\n",
    "print(data.shape)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we were making our CSV file from a pandas dataframe we can use `to_csv`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_df = pd.DataFrame(np.random.random((150, 5)), columns=['A', 'B', 'C', 'D', 'E'])\n",
    "my_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_df.to_csv('randoms.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pandas also provides options for memory mapping the text file to reduce I/O\n",
    "overhead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_df = pd.read_csv('randoms.csv', memory_map=True)\n",
    "my_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performance Tip: Chunking and Iterating\n",
    "\n",
    "So far we've been loading all data at once or depending on memory mapping to\n",
    "reduce the amount of data that was loaded at any one time. Another possible\n",
    "improvement can come from loading chunks of data at a time. This is similar\n",
    "to what we did with the original plain text file iterating over the lines of\n",
    "the file one at a time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By chunking and iterating over the data we can process files that would not\n",
    "fit in memory otherwise. For more info on chunking, see the\n",
    "[pandas documentation](https://pandas.pydata.org/pandas-docs/stable/user_guide/io.html#iterating-through-files-chunk-by-chunk)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reader = pd.read_csv('randoms.csv', iterator=True, chunksize=4)\n",
    "for idx, chunk_df in enumerate(reader):\n",
    "    print(chunk_df)\n",
    "    if idx >= 2: break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This reduces the amount of memory being used at any one time while we work with a single chunk.\n",
    "\n",
    "We can change how big of a chunk we get by calling `get_chunk` with the number of\n",
    "rows to put in the DataFrame returned. Note how we are continuing our reading from the above cells since the reader is an iterator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reader.get_chunk(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CSV - Pros/Cons\n",
    "\n",
    "<details>\n",
    "    <summary>Pros</summary>\n",
    "\n",
    "* Human readable\n",
    "* Can be read in Microsoft Excel (non-programmer collaboration)\n",
    "* Lazy/iterable loading possible\n",
    "* Row-based operations are relatively fast\n",
    "\n",
    "</details>\n",
    "    \n",
    "<details>\n",
    "    <summary>Cons</summary>\n",
    "\n",
    "* Slow to read/write\n",
    "* Wasted disk space\n",
    "* Require reading all columns to get value for a single column\n",
    "* Column-based operations are slow (see below)\n",
    "\n",
    "</details>\n",
    "\n",
    "An unfortunate storage layout for a 3 column CSV (8-bit ASCII characters):\n",
    "\n",
    "<table>\n",
    "    <tr>\n",
    "        <th>Byte Offset</th>\n",
    "        <th>Value</th>\n",
    "        <th>Column</th>\n",
    "    </tr>\n",
    "    <tr><td>0</td><td>0.4433,</td><td>1</td></tr>\n",
    "    <tr><td>7</td><td>0.623,</td><td>2</td></tr>\n",
    "    <tr><td>13</td><td>setosa\\n</td><td>3</td></tr>\n",
    "    <tr><td>20</td><td>0.8866,</td><td>1</td></tr>\n",
    "    <tr><td>27</td><td>0.31,</td><td>2</td></tr>\n",
    "    <tr><td>32</td><td>virginica\\n</td><td>3</td></tr>\n",
    "    <tr><td>42</td><td>0.6644,</td><td>1</td></tr>\n",
    "    </table>\n",
    "\n",
    "The above table shows how a CSV file may be stored on disk. If we don't force all floating point numbers to have the same number of digits or string fields are not all the same size then we can't be sure how to quickly get all values for a single column (by calculating offsets) without reading every value for every column. This makes doing column-based calculations, like the average of an entire field, very slow and wasteful. If we are doing a calculation that requires all values in a single row, then this structure is fairly convenient; we can parse one row at a time efficiently.\n",
    "\n",
    "### Further Reading\n",
    "\n",
    "* [Dask DataFrames](https://docs.dask.org/en/latest/dataframe.html) for distributed and lazy pandas dataframes\n",
    "\n",
    "## Review 1\n",
    "\n",
    "* Store as text: Human-readable but slow\n",
    "* Store as binary: Easily indexable and fast to read/write\n",
    "* Memory maps: Better I/O in most cases\n",
    "* Chunking: Good when working on one chunk at a time\n",
    "* Flat binary: Simple, quick solution, multiple dimensions, no structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Parquet\n",
    "\n",
    "Parquet is another tabular data format and can be thought of as the high\n",
    "performance binary version of the CSV file.\n",
    "Analytics workflows typically don't need to read every column from a tabular\n",
    "format and storing data in something like a CSV file can be very wasteful.\n",
    "The first difference in what Parquet brings to the table (pun intended) is storing data by\n",
    "column (image right) instead of by row (image left).\n",
    "\n",
    "<img src=\"https://www.kdnuggets.com/wp-content/uploads/dremio-table-1.jpg\">\n",
    "<center><sub>Credit: https://www.kdnuggets.com/2017/02/apache-arrow-parquet-columnar-data.html</sub></center>\n",
    "<br>\n",
    "\n",
    "If we keep the columns together it is easier to access one column more easily\n",
    "and more quickly than data stored row by row.\n",
    "\n",
    "## Performance Tip: Spatial and Temporal Locality\n",
    "\n",
    "Modern CPUs have multiple levels of caching. The closer a cache level is to the CPU (the computation) the smaller it is and less it can store at any one time. These closer cache levels are also much faster to get data from. Conversely, caches that are further from the CPU are larger and slower to access. The diagram shows the various caching levels common in modern computers where L1 caches are the smallest and fastest, then L2, L3, main RAM memory, and finally the hard disk; the largest and slowest storage on the local machine. \n",
    "\n",
    "<img width=\"350px\" src=\"https://softwarerajivprab.files.wordpress.com/2019/07/cache.png\" alt=\"CPU L1 L2 L3 cache diagram\">\n",
    "<center><sub>Credit: https://software.rajivprab.com/2018/04/29/myths-programmers-believe-about-cpu-caches/</sub></center><br>\n",
    "\n",
    "If we want the best performance out of the file format we are using then we want to do as many operations as possible with what is in the L1 cache before replacing it with new data. Otherwise, we could suffer from reloading the same data from slower caches. **Temporal locality** is the idea that if we access a particular memory location, we are very likely to access that same location again in the near future. **Spatial locality** is the idea that if we access a particular memory location, an operation in the near future is probably going to involve the data right next to it. Modern computers will assume this to be true and will predictively cache things to get the best performance. That means our best option is to use the memory the way the computer thinks we are going to use it.\n",
    "\n",
    "Further Reading: [Wikipedia](https://en.wikipedia.org/wiki/Locality_of_reference)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performance Tip: Single Instruction, Multiple Data (SIMD)\n",
    "\n",
    "In addition to the multiple cache levels, modern CPUs can also operate on multiple pieces of data at the same time with a single CPU instruction. These vectorized instructions are referred to as SIMD. The CPU is told to do one operation on many values (ex. add 5 to every value) and can perform it on multiple values at a time, in parallel. \n",
    "\n",
    "Even though SIMD instructions are low-level, NumPy does a lot of the work for us by telling the CPU to use SIMD instructions when possible. However, this usually depends on taking advantage of locality (see above) so that the CPU has all the values it is going to operate on.\n",
    "\n",
    "Further Reading: [Wikipedia](https://en.wikipedia.org/wiki/SIMD)\n",
    "\n",
    "<br>\n",
    "\n",
    "Parquet's design for how data is stored on disk and operated on in-memory tries to take advantage of these concepts. Let's run through a basic parquet example to see how we can read and write a parquet file and how we can control some aspects of these complex topics. We'll start by generating some random data to play with. We'll re-assign the `'B'` column so not all of our data is random (to show the results of compression better)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "data = pd.DataFrame(np.random.random((100_000, 5)), columns=['A', 'B', 'C', 'D', 'E'])\n",
    "data['B'] = np.repeat(np.arange(100), 1000)\n",
    "print(data.shape)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's write this data to a local parquet file using the `fastparquet` library (there are others):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fastparquet as fp\n",
    "fp.write('randoms_default.parq', data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls -lh randoms_default.parq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parq_file = fp.ParquetFile('randoms_default.parq')\n",
    "parq_file.to_pandas().head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default fastparquet makes a single file with no compression and only stores entire columns together (no row-groups). Although these defaults may not provide the best performance for a particular use case, the format itself provides us some nice conviences. For example, since the data is stored by column we can quickly and efficiently load a limited set of the columns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data2 = parq_file.to_pandas(['C', 'D'])\n",
    "data2.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take advantage of some of the other features of parquet by writing a new file from our original iris data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp.write('randoms_gzip.parq', data, compression='GZIP')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls -lh randoms*.parq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that using compression can save us a lot of space in our simple data file.\n",
    "\n",
    "So far we've been storing entire columns contiguously, one large block of data. It is also possible to define \"row groups\" to better control how our data is organized on disk. Row groups are grouping a certain number of rows together, but still by column. Let's look at what a parquet file with row groups might look like:\n",
    "\n",
    "<img width=\"400px\" src=\"https://www.dremio.com/img/blog/parquet_block1.png\" alt=\"Parquet basic row group diagram\">\n",
    "<center><sub>Credit: https://www.dremio.com/tuning-parquet/</sub></center>\n",
    "\n",
    "By defining row groups (the orange groups in the above image), we allow ourselves better performance when working on row-based operations while still getting the advantages of storing data by column. We can define our row group sizes using `row_group_offsets` keyword argument in `fastparquet`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp.write('randoms_row_group.parq', data, compression='GZIP', row_group_offsets=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls -lh randoms*.parq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pandas has its own methods for reading parquet files (using the PyArrow library underneath):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "%timeit pd.read_parquet('randoms_row_group.parq').head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And Dask can also load parquet files (using fastparquet underneath). By using Dask we can load the data lazily and in parallel with multiple worker threads or processes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dask.dataframe as dd\n",
    "%timeit dd.read_parquet('randoms_row_group.parq').head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performance Tip: Block sizes\n",
    "\n",
    "Each storage device (cache, RAM, disk) will typically have a default size for the single \"chunks\" or \"units\" of data that they will operate on or store in one contiguous location or provide to another storage element. This is typically referred to as the \"block size\" and can have different names depending on what is being talked about and by whom. To get the best performance we want to try to operate on and store our data in a way that is compatible with the block size of our storage device. If we want to work with 5KB chunks, but our storage device operates on 4KB blocks, we will need to load two 4KB blocks for every 5KB chunk we want to work with. By aligning your workflow with the block size of your storage you can avoid unnecessary I/O operations and delays. However, this isn't always something you have control over or knowlege of.\n",
    "\n",
    "For example, with parquet we want to store our row groups with sizes that take advantage of the block size of the storage while getting the most out of the parquet format. See [Tuning Parquet](https://www.dremio.com/tuning-parquet/) for more details.\n",
    "\n",
    "<img src=\"https://www.dremio.com/img/blog/parquet_block2.png\" alt=\"Parquet row group block sizes\">\n",
    "<center><sub>Credit: https://www.dremio.com/tuning-parquet/</sub></center>\n",
    "\n",
    "<br>\n",
    "\n",
    "The last important parquet feature we should talk about is the ability to store a single parquet dataset as a directory of files instead of one single file. This structure can allow better performance when storing data on cloud file systems or other high performance storage solutions (ex. Apache Hive). We can control this in `fastparquet` by using the `file_scheme` keyword argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp.write('randoms_row_group_hive.parq', data, compression='GZIP', row_group_offsets=1000, file_scheme='hive')\n",
    "!ls -lh randoms_row_group_hive.parq/ | head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parquet - Pros/Cons\n",
    "\n",
    "<details>\n",
    "    <summary>Pros</summary>\n",
    "\n",
    "* Column-based selection very fast\n",
    "* Multiple compression algorithms available\n",
    "* Row groups allow for more control over storage\n",
    "* Cloud storage friendly\n",
    "* Can be stored as single file or directory of files\n",
    "\n",
    "</details>\n",
    "    \n",
    "<details>\n",
    "    <summary>Cons</summary>\n",
    "\n",
    "* Complex (lots of options and choices)\n",
    "* Limited to tabular data (mostly)\n",
    "\n",
    "</details>\n",
    "<br>\n",
    "Further Reading:\n",
    "\n",
    "* [Tuning Parquet](https://www.dremio.com/tuning-parquet/)\n",
    "* [fastparquet Documentation](https://fastparquet.readthedocs.io/en/latest/index.html)\n",
    "* [fastparquet Usage Notes](https://fastparquet.readthedocs.io/en/latest/details.html)\n",
    "* [Apache Arrow](https://arrow.apache.org/)\n",
    "* [Pandas - Scaling to large datasets](https://dev.pandas.io/docs/user_guide/scale.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## HDF5\n",
    "\n",
    "So far we've been dealing with tabular data where our entire dataset can be represented by a 2D array (rows, columns). Sometimes your data may require a more complex hierarchical structure. One good solution for these more complex structures is the Hierarchical Data Format (HDF); specifically the HDF5 file format. HDF5 is a self-describing format which means that in addition to the raw data you can store all related metadata (when it was recorded, where it was recorded, etc).\n",
    "\n",
    "Pandas and other python libraries do provide utility functions for writing and reading HDF5 files (.h5 extension), but we will be using the `h5py` library here to step through the structure of the files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The HDF5 data model is made up of a few key components or ideas:\n",
    "\n",
    "* File: The single object stored on disk\n",
    "* Datasets: Multidimensional data array with attributes and other metadata\n",
    "* Groups: A collection of objects (datasets or groups)\n",
    "* Attributes: Key-value store of metadata on groups or datasets\n",
    "\n",
    "<img width=\"700px\" src=\"https://www.neonscience.org/sites/default/files/images/HDF5/hdf5_structure3.jpg\" alt=\"Complex HDF5 layout\">\n",
    "<center><sub>Credit: https://www.neonscience.org/about-hdf5</sub></center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "h = h5py.File('test.h5', 'w')\n",
    "h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HDF5 files have an implicit \"root\" group where groups or datasets can be added or attributes can be set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h['/']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.attrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dict(h.attrs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "h.attrs['start_time'] = datetime.utcnow().isoformat(\" \")\n",
    "h.attrs['source'] = \"Dave's Notebook\"\n",
    "h.attrs['revision'] = 1\n",
    "h.attrs['inputs'] = ['a.txt', 'b.txt', 'c.text']\n",
    "dict(h.attrs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can store most builtin Python data types in our attributes to store various pieces of metadata about our file. These are typically referred to as \"global\" attributes.\n",
    "\n",
    "We can also add data to our file as Datasets. A Dataset can be created with the data specified:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds1 = h.create_dataset('dataset1', data=np.random.random((10, 100, 100)))\n",
    "ds1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds1.attrs['valid_range'] = [0, 1]\n",
    "ds1.attrs['missing_value'] = np.nan\n",
    "dict(ds1.attrs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A dataset can also be created and have data filled in later. There are a lot of options available when creating a dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.create_dataset?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds2 = h.create_dataset('dataset2', shape=(5, 100), dtype=np.uint8,\n",
    "                       compression='gzip', chunks=(2, 100), fillvalue=255)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds2[2:4] = np.arange(100)\n",
    "ds2[0] = np.arange(100, 200)\n",
    "ds2[-1] = np.arange(150, 250)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whether we are reading or writing we can use slicing syntax to access parts of the data array (only loading/writing what's needed) or the entire thing.\n",
    "\n",
    "If we look at the data we wrote to the `\"dataset2\"` dataset, notice how the second row of data uses the `fillvalue` we set above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds2[:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can create groups to describe more complex structures of our data (groups with groups with datasets, etc):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1 = h.create_group('geolocation')\n",
    "g1.attrs['my_key'] = 'my_val'\n",
    "# sub-group of g1\n",
    "g12 = g1.create_group('geocentric')\n",
    "# another group off of the root\n",
    "g2 = h.create_group('observed_data')\n",
    "# dataset as part of a group\n",
    "obs_a = g2.create_dataset('observation_a', data=np.arange(10.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The HDF5 C library comes with the `h5dump` command line tool for investigating an HDF5 file in different ways. Here we use the `-n` flag to list the contents of the file without their attributes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!h5dump -n test.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HDF5 - Pros/Cons\n",
    "\n",
    "<details>\n",
    "    <summary>Pros</summary>\n",
    "\n",
    "* Complex structures\n",
    "* Per-dataset compression\n",
    "* Multidimensional arrays\n",
    "* Metadata\n",
    "\n",
    "</details>\n",
    "    \n",
    "<details>\n",
    "    <summary>Cons</summary>\n",
    "\n",
    "* Complex (lots of options and choices)\n",
    "* Limited parallel/thread-safe operations\n",
    "* Not cloud storage friendly\n",
    "\n",
    "</details>\n",
    "<br>\n",
    "\n",
    "HDF5 supports a lot of complex features including references HDF5 objects from external files or using references to point to other locations in the same file. These features are left as an exercise for the reader to explore.\n",
    "\n",
    "Further Reading:\n",
    "\n",
    "* [HDF5 Data Model](https://support.hdfgroup.org/HDF5/doc1.6/UG/03_Model.html)\n",
    "* [h5py Documentation](http://docs.h5py.org/en/stable/index.html)\n",
    "* [HDF5 Reference Objects](http://docs.h5py.org/en/stable/refs.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## NetCDF4\n",
    "\n",
    "NetCDF4 is another format that supports storing multidimensional arrays like HDF5, but provides a simpler/flatter structure to ease use of the files. Under the hood NetCDF4 is actually built on top of the HDF5 format. As such, almost all of the features provided by HDF5 are also available in NetCDF4 (chunking, per-dataset compression, etc). NetCDF4 is a very popular and common format for storing scientific data. When mixed with metadata standards like those defined by the [Climate and Forecast (CF) group](http://cfconventions.org/) files become much easier to distribute and use between different applications.\n",
    "\n",
    "The NetCDF4 C library not only allows you to read and write NetCDF4 files, but can act as a generic reading library for formats that fulfill the NetCDF4 Data Model. This means that the NetCDF4 C library can read NetCDF4, NetCDF3, and HDF4 files. In future versions it will also be able to read the Zarr format (see next section). The NetCDF4 data model is very similar to HDF5, but with a few key differences. The first major difference is that HDF5 Datasets are called \"Variables\". The second major addition is the more formal definition of labeled Dimensions.\n",
    "\n",
    "<img src=\"https://www.unidata.ucar.edu/software/netcdf/docs/nc4-model.png\" alt=\"NetCDF4 Data Model\">\n",
    "<center><sub>Credit: https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_data_model.html</sub></center>\n",
    "\n",
    "We'll start exploring the NetCDF4 data model by creating a simple NetCDF4 file using the `netcdf4-python` library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import netCDF4\n",
    "\n",
    "nc = netCDF4.Dataset('test.nc', 'w')\n",
    "nc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll add some data, but before we do that we need to define the dimensions for our variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_dim = nc.createDimension('y')\n",
    "x_dim = nc.createDimension('x', 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var1 = nc.createVariable('var1', np.uint16, dimensions=('y', 'x'))\n",
    "var1[:] = np.random.randint(1, 65535, size=(200, 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also define variables with the same name as a dimension. These are known as \"coordinate\" variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_var = nc.createVariable('x', np.float32, dimensions=('x',))\n",
    "x_var[:] = np.linspace(1000.0, 1500.0, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like HDF5 we can assign metadata attributes to our variables or to the global attributes of the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nc.created = 'now'\n",
    "var1.valid_range = [0, 65535]\n",
    "var1.standard_name = 'toa_brightness_temperature'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we defined the `y` dimension we didn't specify a size. In NetCDF4 this means the dimension is 'unlimited' and can grow as needed. We can check the size of a dimension by looking at `len(dim)`. Let's add more data to the `var1` variable and see how the dimensions change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(y_dim))\n",
    "var1[205] = np.arange(100)\n",
    "print(len(y_dim))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nc['var1']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nc.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the easiest ways to work with NetCDF files in python is through the Xarray library by using it's `open_dataset` and `to_netcdf` functions. Xarray can take advantage of the dask library for lazy loading our data and computing things in parallel. By specifying the `chunks` keyword argument we can tell Xarray to use dask.\n",
    "\n",
    "<img width=\"700px\" src=\"http://xarray.pydata.org/en/stable/_images/dataset-diagram.png\">\n",
    "<center><sub>Credit: http://xarray.pydata.org/en/stable/data-structures.html</sub></center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "xnc = xr.open_dataset('test.nc', chunks={'y': 50, 'x': 100})\n",
    "xnc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(xnc.dims)\n",
    "print(xnc.attrs)\n",
    "print(xnc['var1'].dims)\n",
    "print(xnc['var1'].coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xnc['var1']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Xarray helps us when slicing data by applying the same slices to our coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xnc['var1'][:50, 50:75]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xnc.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to `h5dump` the NetCDF4 library comes with a `ncdump` utility:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ncdump -h test.nc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NetCDF4 - Pros/Cons\n",
    "\n",
    "<details>\n",
    "    <summary>Pros</summary>\n",
    "\n",
    "* Simpler HDF5 structures\n",
    "* Per-variable compression and chunking\n",
    "* Multidimensional arrays\n",
    "* Named dimensions and coordinate variables\n",
    "* Data model compatible with Xarray (limited support for groups)\n",
    "* Metadata\n",
    "\n",
    "</details>\n",
    "    \n",
    "<details>\n",
    "    <summary>Cons</summary>\n",
    "\n",
    "* Limited parallel/thread-safe operations\n",
    "* Not cloud storage friendly\n",
    "\n",
    "</details>\n",
    "<br>\n",
    "\n",
    "Further Reading:\n",
    "\n",
    "* [xarray.open_mfdatasets](http://xarray.pydata.org/en/stable/generated/xarray.open_mfdataset.html) - Load multiple files as one Dataset\n",
    "* [OpenDAP](https://en.wikipedia.org/wiki/OPeNDAP) - Read remote NetCDF files without loading the entire file\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Zarr\n",
    "\n",
    "One of the major downsides to HDF5 and NetCDF4 is that they do not allow for parallel processing or do not work well with cloud storage solutions. To address these limitations and others the Zarr format was developed. Zarr is developed heavily by the scientific Python community including the developers of dask and xarray. In short, Zarr is like a NetCDF4 file if it was stored as a directory of files instead of one single file.\n",
    "\n",
    "Since Zarr uses the same data model as NetCDF/Xarray, we can skip directly to what the format looks like and how it is used. We'll start by using the NetCDF file we made before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "xnc = xr.open_dataset('test.nc', chunks={'y': 50, 'x': 100})\n",
    "xnc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xnc.to_zarr('test.zarr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tree -a test.zarr/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above we can see that `to_zarr` defaulted to using the chunking of our dask arrays when writing to zarr storage. Each chunk shows up as its own file with metadata (both user metadata and array information) as separate files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat test.zarr/var1/.zarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat test.zarr/var1/.zattrs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To read zarr files we can either use the `zarr` python library or use Xarray again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zds = xr.open_zarr('test.zarr')\n",
    "zds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Zarr - Pros/Cons\n",
    "\n",
    "<details>\n",
    "    <summary>Pros</summary>\n",
    "\n",
    "* Designed for efficient storage and use\n",
    "* Per-chunk compression\n",
    "* Multiple storage options (zip file, database, etc)\n",
    "* Multidimensional arrays\n",
    "* Named dimensions and coordinate variables\n",
    "* Data model compatible with Xarray (limited support for groups)\n",
    "* Metadata\n",
    "* Future support in NetCDF4 C library\n",
    "\n",
    "</details>\n",
    "    \n",
    "<details>\n",
    "    <summary>Cons</summary>\n",
    "\n",
    "* Relatively new\n",
    "* Library not available in every popular programming language\n",
    "\n",
    "</details>\n",
    "<br>\n",
    "\n",
    "Further Reading:\n",
    "\n",
    "* [zarr Documentation](https://zarr.readthedocs.io/en/stable/)\n",
    "* [Builtin zarr Storages](https://zarr.readthedocs.io/en/stable/api/storage.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "* Text is slow, but might be needed if humans need to read it\n",
    "* Text is the easiest format to share with other languages/programs (excel)\n",
    "* Compression is good for space saving but not read/write speeds\n",
    "* Chunking data is good\n",
    "* Parquet and Zarr are cloud friendly\n",
    "* File formats are complex\n",
    "* There are so many file formats\n",
    "\n",
    "This Notebook:\n",
    "https://github.com/djhoese/data-file-formats\n",
    "\n",
    "## Honarable Mentions\n",
    "\n",
    "These files are no less useful than the above formats, but we have to draw the line somewhere.\n",
    "\n",
    "* Cloud Optimized GeoTIFFs (COG)\n",
    "* Apache Avro (row-based)\n",
    "* Apache Arrow (in-memory)\n",
    "* ORC\n",
    "* JSON\n",
    "* XML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rm -r *.zarr *.parq *.csv *.dat *.txt *.nc *.h5 *.npy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.3"
  },
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}