01_Basic_Data_Analysis_and_Visualization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesson 1: Basic Data Analysis and Visualization \n",
    "\n",
    "Author: Rebekah Esmaili (rebekah.esmaili@gmail.com)\n",
    " \n",
    "---\n",
    "\n",
    "\n",
    "## Why Python?\n",
    "\n",
    "Pros\n",
    "\n",
    "* General-purpose, cross-platform\n",
    "* Free and open source\n",
    "* Reasonably easy to learn\n",
    "* Expressive and succinct code, forces good style\n",
    "* Being interpreted and dynamically typed makes it great for data analysis\n",
    "* Robust ecosystem of scientific libraries, including powerful statistical and visualization packages\n",
    "* Large community of scientific users and large existing codebases\n",
    "* Major investment into Python ecosystem by Earth science research agencies, including NASA, NCAR, UK Met Office, and Lamont-Doherty Earth Observatory. See Pangeo.\n",
    "* Reads Earth science data formats like HDF, NetCDF, GRIB\n",
    "\n",
    "Cons\n",
    "\n",
    "* Performance penalties for interpreted languages, although many libraries are wrappers for compiled languages. Avoid large loops in favor of matrix/vector operations when possible.\n",
    "* Multithreading is limited due to the Global Interpreter Lock, but other parallelism is available\n",
    "* See Julia for a modern scientific language which is trying to overcome these challenges\n",
    "\n",
    "Why you should not use Python 2.7:\n",
    "\n",
    "* Python 2 reached it's \"end of life\" as of January 2020\n",
    "* No more updates or bugfixes\n",
    "* No further official support\n",
    "* Subtle differences: https://www.geeksforgeeks.org/important-differences-between-python-2-x-and-python-3-x-with-examples/\n",
    "\n",
    "---\n",
    "\n",
    "## Lesson Objectives\n",
    "\n",
    "* You will learn to:\n",
    "    * Import relevant packages for scientific programming\n",
    "    * Read ascii data\n",
    "    * Basic plotting and visualization\n",
    "   \n",
    "---\n",
    "\n",
    "## Basic Python Syntax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most basic Python command is to write words to the screen. In jupyter notebooks, the result will appear below the line of code. To run the above command in Jupyter notebook, highlight the cell and either chick the run button (►) or press the **Shift** and **Enter** keys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a comment, python will not run this!\n",
    "print(\"Hello Earth\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Python, variables are dynamically allocated, which means that you do not need to declare the type or size prior to storing data in them. Instead, Python will automatically guess the variable type based on the content of what you are assigning:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_int = 8\n",
    "var_float = 15.0\n",
    "var_scifloat = 4e8\n",
    "var_complex = complex(4, 2)\n",
    "var_greetings = 'Hello Earth'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python has many built in functions, the syntax is usually:\n",
    "\n",
    "```\n",
    "function_name(inputs)\n",
    "```\n",
    "You have already used two functions: *print()* and *complex()*. Another useful function is *type()*, will tell us if the variable is an integer, a float, a complex number, or a string. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " type(var_int),  type(var_float), type(var_scifloat), type(var_complex), type(var_greetings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python has the following built-in operators:\n",
    "\n",
    "* Addition, subtraction, multiplication, division: +, -, *, /\n",
    "* Exponential, integer division, modulus: \\**, //, %\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "2+2.0, var_int**2, var_float//var_int, var_float%var_int"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "**Exercise 1:** Learning to use notebooks\n",
    "\n",
    "1. Launch Jupyter Notebook and create a new notebook\n",
    "2. Rename the notebook\n",
    "3. Create a new cell and use *type()* to see if the following are floats and integers:\n",
    "    * 2+2\n",
    "    * 2\\*2.0\n",
    "    * var_float/var_int\n",
    "---\n",
    "**Solution:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with lists\n",
    "\n",
    "Lists are useful for storing scientific data. Lists are made using square brackets. They can hold any data type (integers, floats, and strings) and even mixtures of the two."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numbers_list = [4, 8, 15, 16, 23]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can access elements of the list using the index. Python is zero based, so index 0 retrieves the first element."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numbers_list[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "New items can also be appended to the list using the append function, which has the syntax:\n",
    "\n",
    "```\n",
    "variable.function(element(s))\n",
    "```\n",
    "The list will be updated *in-place*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numbers_list.append(42)\n",
    "print(numbers_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perhaps we want to calculate the sum of the values in two lists. However, we cannot use the *+* like we did with single values. For list objects, the + will *combine* lists."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numbers_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To perform mathematical operations, you can convert the above list to an array using the NumPy package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing Packages\n",
    "\n",
    "Packages are collection of modules, which help simplify common tasks. [NumPy](https://numpy.org/) is useful for mathematical operations and array manipulation.\n",
    "\n",
    "\n",
    "* Provides a high-performance multidimensional array object and tools for working with these arrays.\n",
    "* Fundamental package for scientific computing with Python.\n",
    "* Included with with the Anaconda package manager.\n",
    "* For more examples than presented below, please refer [the NumPy Quick Start](https://numpy.org/devdocs/user/quickstart.html)\n",
    "\n",
    "The basic syntax for calling packages is to type the import \\[package name\\]. However, some packages have long names, so you can use import \\[package name\\] as \\[alias\\]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you do not see any error after running the line above, then the package was successfully imported.\n",
    "### Working with arrays\n",
    "\n",
    "I can use NumPy’s array constructor *np.array()* to convert our list to a NumPy array and perform the matrix multiplication. For example, I can double each element of the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numbers_array = np.array(numbers_list)\n",
    "numbers_array*2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another difference between arrays and lists is that lists are only one-dimensional. NumPy can be any number of dimensions. For example, I can change the dimensions of the data using the *reshape()* function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numbers_array_2d = numbers_array.reshape(3,2)\n",
    "numbers_array_2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numbers_array_2d.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The original numbers_array had 6 elements in a 1D array, the new array is 2D with 3 rows and 2 columns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading ASCII data\n",
    "\n",
    "The Pandas package has a useful function for reading text/ascii data called *read_csv()*. The function name is somewhat a misnomer, as *read_csv* will read any delimited data using the *delim=* keyword argument. Below, you will import the [Pandas](https://pandas.pydata.org/) package and we will read in a dataset. Note that the path below is relative to the current notebook and you may have to change the code if you are running locally on your computer:\n",
    "\n",
    "```\n",
    "data/VIIRSNDE_global2020258.v1.0.txt\n",
    "```\n",
    "\n",
    "We will look at the Visible Infrared Imaging Radiometer Suite (VIIRS) Active Fire product, a product that classifies if a pixel contains fire with various confidence levels. More information can be found at https://www.ospo.noaa.gov/Products/land/fire.html. We will examine the data on Sept 15, 2020 (day of year 258)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default seperator is a comma (,), however my data also contains space. I use the \"\\s*\" to indicate space following the comma should be ignored. The engine=\"python\" keyword ensures that this will work across different operating systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = \"data/VIIRSNDE_global2020258.v1.0.txt\"\n",
    "fires = pd.read_csv(fname, sep=',\\s*', engine='python')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can inspect the contents within the notebook using the *head()* function, which will return the first five rows of the dataset. Pandas automatically stores data in structures called *DataFrames*. DataFrames are two dimensional (rows and columns) and resemble a spreadsheet. The leftmost column is the row index and is not part of the *fires* dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fires.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can access individual columns of data using the column name. For example, below you can extract the pixel brightness temperature (brt):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fires[\"brt_t13(K)\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "**Exercise 2:** Import an ascii file\n",
    "\n",
    "1. Import the dataset \"20200901_20200930_Monterey.lev15.csv\" and save it to a variable called *aeronet*.\n",
    "2. Print the first few lines using *.head()*\n",
    "3. Find a column that doesn't have only missing values (-999) and (challenge!) calculate the mean using the following syntax *variable\\[\"column\"\\].mean()*\n",
    "---\n",
    "**Solution:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Working with masks and masked arrays\n",
    "\n",
    "When working with data, sometimes there are numbers I want to remove. For instance, I may want to work with data below a certain threshold. You can subset the data using identity operations:\n",
    "\n",
    "* less than: <\n",
    "* less than or equal to: <=\n",
    "* greater than: >\n",
    "* greater than or equal to: >=\n",
    "* equals: ==\n",
    "* not equals: !=\n",
    "\n",
    "Their use will return either a True or False statement. For the *fires* dataset, you can find which elements of the array that meet some condition, such as only examining larger fires that have a Fire Radiative Power (FRP) above 50 MW:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "masked_nums = (fires['frp(MW)'] > 50)\n",
    "print(masked_nums)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes you may want to filter by two conditions. For example, insteading of filtering the FRP data, you may only want to examine values within a latitude and longitude domain. In Python, I can combine multiple conditions using and (&) and or (|) statements. Below, I extract the data in 5&deg;x5&deg; box arond Monterey, California:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "masked_nums = (fires['Lat'] > 35.0) & (fires['Lat'] < 40.0) & (fires['Lon'] > -125.0) & (fires['Lon'] < -120.0)\n",
    "print(masked_nums)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above mask can be used in place of an index. Below, you can create a new variable that takes the FRP using the *fires\\['frp(MW)'\\]* variable and subsets it with the array of *masked_nums*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "monterey_fires = fires['frp(MW)'][masked_nums]\n",
    "print(monterey_fires)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this new variable, you can compute the average in this region and compare them to the global average for that day:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "monterey_fires.mean(), fires['frp(MW)'].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use the size command to compare the dimensions of original array and the one that filtered out values that were outside of our latitude and longitude bounds. You will notice that these two arrays have different sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fires['frp(MW)'].size, monterey_fires.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are cases where you will want to preserve the size and shape of the original array. For these situations, you can utilize the NumPy *masked array* module. The syntax is *np.ma.array()*, and you will add a keyword argument *mask=*, which is set to the inverse (~) of the *mask_nums*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "monterey_fires_ma = np.ma.array(fires['frp(MW)'], mask=~masked_nums, fill_value=-999)\n",
    "monterey_fires_ma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, you can calculate the mean values and confirm that they are the same as the previous example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "monterey_fires_ma.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, the key difference will be the size, which retains the shape of the unmasked data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "monterey_fires_ma.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**Exercise 3:** Filtering data\n",
    "\n",
    "Using the dataset imported in the previous example (*aeronet*):\n",
    "    \n",
    "1. Create a mask that filters the \"AOD_870nm\" column to only include values that are above 0.\n",
    "2. Create a new variables, *day_of_year*, with the mask applied to aeronet\\[\"Day_of_Year(Fraction)\"\\].\n",
    "3. Create a new variables, *aod_870*, with the mask applied to aeronet\\[\"AOD_870nm\"\\].\n",
    "4. Compare the mean value of *aeronet\\[\"AOD_870nm\"\\]* to *aod_870*.\n",
    "    \n",
    "---\n",
    "**Solution**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basic figures and plots\n",
    "\n",
    "Python has several packages to create visuals for remote sensing data, either in the form of imagery or plots of relevant analysis. Of these, the most widely used and oldest packages is [Matplotlib](https://matplotlib.org/). Matplotlib plots are highly customizable and has additional toolkits that can enhance functionality, such as creating maps using the [Cartopy](https://scitools.org.uk/cartopy/docs/latest/) package, which I will describe more in the next session.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose you want to learn what the global distribution of fire radiative power is. From inspecting the frp(MW) column earlier, these values extend to many decimal places. Rather than use a continuous scale, I can instead group in the data into 10 MW bins, from 0 to 500 MW:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins10MW = np.arange(0, 500, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I can use these bins to create a histogram. Line by line, the code below will do as follows. Each additional line is layering elements on this empty graphic. The entire block of code must be run at once and not split into multiple cells.  \n",
    "\n",
    "1. *plt.figure()* creates a blank canvas.\n",
    "2. I add the histogram to the figure using *plt.hist()*, which automatically will count the number of rows with fire radiative power in the bins that I defined above in the bins10W variable. I must then pass in the data (fires['frp(MW)']) and the bins (bins10MW) into plt.hist. \n",
    "3. *plt.show()* tells matplotlib the plot is now complete and to render it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=[5,5])\n",
    "plt.hist(fires['frp(MW)'], bins=bins10MW)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, you will remake this plot but add some aesthetic additions, such as labels to the x and y axis using *set_xlabel()* and *set_ylabel()*. Since there are thousands more fires with fire radiative power less than 100 MW than fires with higher values the data are likely lognormal. The plot will be easier to interpret of I rescale the y-axis to a log scale while leaving the x-axis linear.\n",
    "\n",
    "The command *plt.subplot()* will return an axis object to a variable (*ax*). There are three numbers passed in (111), which correspond to rows, columns, and index. In this example, there is one row and one column, and therefore, only one index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "\n",
    "ax = plt.subplot(111)\n",
    "\n",
    "ax.hist(fires['frp(MW)'], bins=bins10MW)\n",
    "\n",
    "ax.set_yscale('log')\n",
    "\n",
    "ax.set_xlabel(\"Fire Radiative Power (MW)\")\n",
    "ax.set_ylabel(\"Counts\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also plot the data in 2-dimensions. For example, each row in *fires* has a latitude and longitude coordinates pair. I will take these two coordinates and plot using *plt.scatter()*. The first argument is the x-coordinate and the second is the y-coordinate (the order matters). \n",
    "\n",
    "There are some command line options *plt.scatter()*:\n",
    "\n",
    "* s: size with respect to the default\n",
    "* c: color, which can be either from a predefined name list or a hexadecimal value\n",
    "* alpha: opacity, where smaller values are transparent.\n",
    "\n",
    "Like in the previous example, I have chosen to label the latitude and longitude axes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.subplot(111)\n",
    "\n",
    "ax.scatter(fires['Lon'], fires['Lat'], s=0.5, c='black', alpha=0.1)\n",
    "\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can almost see the outline of the continents from the data above. In the next session, you will learn how to overlay maps onto your plots."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**Exercise 4:** Create a scatterplot\n",
    "\n",
    "Use the variables *aod_870* and *day_of_year* that you made in Exercise 3 to:\n",
    "\n",
    "1. Create a scatter plot showing the *day_of_year* (x-axis) and *aod_870* (y-axis)\n",
    "2. Add y-axis and x-axis labels using *.set_xlabel()* and *.set_ylabel()*\n",
    "3. Adjust the color and size of the scatterplot\n",
    "---\n",
    "**Solution**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary:\n",
    "\n",
    "You learned:\n",
    "* Very basic built-in Python functions and operations\n",
    "* How to import three packages: numpy, pandas, and matplotlib\n",
    "* Worked with arrays and lists\n",
    "* How to create a simple plot\n",
    "\n",
    "Next lesson:\n",
    "* More advanced plots, such as using maps\n",
    "* Importing scientific datasets, such as netcdf and grib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.11.0 ('notebook_demo')",
   "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.11.0"
  },
  "vscode": {
   "interpreter": {
    "hash": "4589143d4cda0c8671911bd60c16dc1d10ec327722e7574bc882b745b51509b4"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}