{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "950099de-bfc3-4e1d-85a4-0184030e85a8",
   "metadata": {},
   "source": [
    "<table>\n",
    "    <tr>\n",
    "        <td>\n",
    "            <img src=\"https://arm.gov/uploads/TRACER_skyline.jpg\" style=\"height: 200px;\" alt=\"TRACER Campaign\"></img>\n",
    "        </td>\n",
    "        <td>\n",
    "            <h1>Atmospheric Radiation Measurement user facility (ARM)</h1>\n",
    "            <h1>TRacking Aerosol Convection interations ExpeRiment (TRACER)</h1>\n",
    "            <h2>July 16-19, 2022 Dust Event</h2>\n",
    "            Notebook for data exploration of ARM aerosol and lidar data. </br>\n",
    "            Corresponding Author: Adam Theisen (atheisen@anl.gov)\n",
    "        </td>\n",
    "    </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "830e2e64-3205-41cd-8247-2fa3905f6106",
   "metadata": {},
   "source": [
    "## Overview\n",
    "The first notebook in this series was an introduction to some of the features in ACT.  In this notebook, we will be exploring different types of data and bringing it all together to get a view of this dust event.\n",
    "\n",
    "1. Micropulse Lidar (MPL) Data\n",
    "1. Aerodynamic Particle Sizer (APS) Data\n",
    "1. Putting it All Together\n",
    "1. Questions for the User to Explore"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74e672b9-ef14-4f01-be19-a222f5ea1573",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Prerequisites\n",
    "This notebook will rely heavily on Python and the [Atmospheric data Community Toolkit (ACT)](https://github.com/ARM-DOE/ACT).  Don't worry if you don't have experience with either, this notebook will walk you though what you need to know.\n",
    "\n",
    "You will also need an account and token to download data using the ARM Live webservice.  Navigate to the [webservice information page](https://adc.arm.gov/armlive/) and log in to get your token.  Your account username will be your ARM username.\n",
    "\n",
    "| Concepts | Importance | Notes |\n",
    "| --- | --- | --- |\n",
    "| [ACT](https://github.com/ARM-DOE/ACT) | Helpful | |\n",
    "| 3-ACT-Basics.ipynb | Helpful ||\n",
    "\n",
    "- **Time to learn**: 60 Minutes\n",
    "- **System requirements**:\n",
    "    - Python 3.11 or latest\n",
    "    - ACT v2.0.0 or latest\n",
    "    - numpy\n",
    "    - xarray\n",
    "    - matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f89921ee-514a-4e11-859a-4858ed7354ab",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0de3c05a-6e49-4859-850b-32c989f2cc33",
   "metadata": {},
   "source": [
    "## Imports\n",
    "Let's get started with some data!  But first, we need to import some libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5b00201-8c78-4ef5-b02f-57e55021f5d6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import act\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import glob"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3647028a-cf78-436a-be33-6a0148474bf6",
   "metadata": {},
   "source": [
    "## Micropulse Lidar (MPL) Data\n",
    "ARM has many value-added products (VAPs) that apply corrections, additional quality control, retrievals, and more.  These are a great way to get the data if you don't want to do the added processing yourself.  In this case, there is a VAP for the MPL with the datastream hou30smplcmask1zwangM1.c1 that you could utilize.  For this case though, let's download the instrument data and correct it ourselves using ACT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dc01deb-8e97-4e45-9efb-196b37ca7386",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Set your username and token here!\n",
    "username = 'YourUsername'\n",
    "token = 'YourToken'\n",
    "\n",
    "# Set the datastream and start/enddates\n",
    "datastream = 'houmplpolfsM1.b1'\n",
    "startdate = '2022-07-16'\n",
    "enddate = '2022-07-16'\n",
    "\n",
    "# Use ACT to easily download the data.  Watch for the data citation!  Show some support\n",
    "# for ARM's instrument experts and cite their data if you use it in a publication\n",
    "result = act.discovery.download_arm_data(username, token, datastream, startdate, enddate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acdaa696-ebb5-4bbe-8d1e-a4c66f524964",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Let's read in the data using ACT and check out the data\n",
    "ds_mpl = act.io.read_arm_netcdf(result)\n",
    "\n",
    "# Now we can correct the data\n",
    "ds_mpl = act.corrections.correct_mpl(ds_mpl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cd2407a-7f6c-4e86-b677-e1da101f1786",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Let's plot up the data to see what it looks like. But first,\n",
    "# if you look at the variables, you would see the the variable\n",
    "# we're going to plot has range_bins as it's 2nd dimension.\n",
    "# We want it to be height so we have to swap some coordinates around\n",
    "ds_mpl.coords['height'] = ds_mpl.height\n",
    "ds_mpl = ds_mpl.swap_dims({'range_bins': 'height'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "477de074-9465-47eb-a658-decf6768b077",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# It's a lot of data, so it will take some time!\n",
    "display = act.plotting.TimeSeriesDisplay(ds_mpl, figsize=(10, 6))\n",
    "display.plot('signal_return_co_pol', cvd_friendly=True, vmin=-20, vmax=20)\n",
    "\n",
    "# This instrument collects data over 30 km so let's constrain\n",
    "# this to just the lower 5 km\n",
    "display.set_yrng([0, 5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0a1cd76-d762-4d28-b125-ee203f1e146e",
   "metadata": {},
   "source": [
    "### Polarized Lidar\n",
    "The micropulse lidar is a polarized lidar which means that it can be used to inform us about the shape of the particles that the light is interacting with.  This website from the [Leibniz Institute for Troposheric Research](https://www.tropos.de/en/research/projects-infrastructures-technology/technology-at-tropos/remote-sensing/polarization-lidar) has a nice explanation of polarized lidar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73c9bcfd-7330-4126-b50d-50556c8c8532",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Let's look at the cross polarization to get some more information\n",
    "# It's a lot of data, so it will take some time!\n",
    "display = act.plotting.TimeSeriesDisplay(ds_mpl, figsize=(10, 6))\n",
    "display.plot('signal_return_cross_pol', cvd_friendly=True, vmin=-20, vmax=20)\n",
    "\n",
    "# This instrument collects data over 30 km so let's constrain\n",
    "# this to just the lower 5 km\n",
    "display.set_yrng([0, 5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a93b725a-3216-4eca-8ad0-d9e9e31b1372",
   "metadata": {},
   "source": [
    "## Aerosol Particle Size Data\n",
    "Dust particles are generally under 10 µm or 10,000 nm and ARM has a variety of instruments that measure across different ranges. The chart below shows the ranges for each instrument and where we can expect to see the dust particles.\n",
    "\n",
    "<img src=\"./images/aerosol_sizing.png\" width=900>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2991b44-a831-4879-b655-e2d7e5fef28b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Download the data as before\n",
    "datastream = 'houaosapsM1.b1'\n",
    "result = act.discovery.download_arm_data(username, token, datastream, startdate, enddate)\n",
    "\n",
    "# and read it in\n",
    "ds_aps = act.io.arm.read_arm_netcdf(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f2c2c57-09c9-44d6-8098-8a414ba0fb0d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Let's plot up the size distribution and number concentration\n",
    "# These data are taken every second so it will take a little bit\n",
    "display = act.plotting.TimeSeriesDisplay(ds_aps, subplot_shape=(2,), figsize=(10,8))\n",
    "display.plot('total_N_conc', subplot_index=(0,))\n",
    "display.day_night_background(subplot_index=(0,))\n",
    "display.plot('dN_dlogDp', subplot_index=(1,), vmax=40, cvd_friendly=True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc59f297-718c-4157-9acc-a44a4b801a91",
   "metadata": {},
   "source": [
    "## Putting it All Together\n",
    "Let's start pulling these data together into the same plots so we can see what's going on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20da29fc-2009-4c7d-9490-301d99eb2f20",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Let's start fresh and get all the data and maybe expand the dates\n",
    "startdate = '2022-07-10'\n",
    "enddate = '2022-07-20'\n",
    "\n",
    "# MPL - Let's try the VAP this time\n",
    "datastream = 'hou30smplcmask1zwangM1.c1'\n",
    "result = act.discovery.download_arm_data(username, token, datastream, startdate, enddate)\n",
    "ds_mpl = act.io.arm.read_arm_netcdf(result)\n",
    "\n",
    "# APS\n",
    "datastream = 'houaosapsM1.b1'\n",
    "result = act.discovery.download_arm_data(username, token, datastream, startdate, enddate)\n",
    "ds_aps = act.io.arm.read_arm_netcdf(result)\n",
    "\n",
    "# SMPS\n",
    "datastream = 'houaossmpsM1.b1'\n",
    "result = act.discovery.download_arm_data(username, token, datastream, startdate, enddate)\n",
    "ds_smps = act.io.arm.read_arm_netcdf(result)\n",
    "\n",
    "#ACSM\n",
    "datastream = 'houaosacsmM1.b2'\n",
    "result = act.discovery.download_arm_data(username, token, datastream, startdate, enddate)\n",
    "ds_acsm = act.io.arm.read_arm_netcdf(result)\n",
    "\n",
    "#PSAP\n",
    "datastream = 'houaospsap3w1mM1.b1'\n",
    "result = act.discovery.download_arm_data(username, token, datastream, startdate, enddate)\n",
    "ds_psap = act.io.arm.read_arm_netcdf(result)\n",
    "\n",
    "#SP2\n",
    "datastream = 'houaossp2bc60sM1.b1'\n",
    "result = act.discovery.download_arm_data(username, token, datastream, startdate, enddate)\n",
    "ds_sp2 = act.io.arm.read_arm_netcdf(result)\n",
    "\n",
    "# MET\n",
    "datastream = 'houmetM1.b1'\n",
    "result = act.discovery.download_arm_data(username, token, datastream, startdate, enddate)\n",
    "ds_met = act.io.arm.read_arm_netcdf(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b2770c8-f40f-4847-b426-02d69bd9529b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Before we proceed to plotting, let's reduce the MPL data down a little bit\n",
    "# This will remove all data where heights are greater than 5\n",
    "ds_mpl = ds_mpl.where(ds_mpl.height <= 5, drop=True)\n",
    "\n",
    "# This will resample to 1 minute\n",
    "ds_mpl = ds_mpl.resample(time='1min').nearest()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db93b1f6-a836-4c28-b7a9-87053d5f9963",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Let's not forget about QCing the data!\n",
    "# We can remove all the bad data from each aerosol dataset\n",
    "ds_aps.clean.cleanup()\n",
    "ds_aps = act.qc.arm.add_dqr_to_qc(ds_aps)\n",
    "ds_aps.qcfilter.datafilter(rm_assessments=['Bad'], del_qc_var=False)\n",
    "\n",
    "ds_smps.clean.cleanup()\n",
    "ds_smps = act.qc.arm.add_dqr_to_qc(ds_smps)\n",
    "ds_smps.qcfilter.datafilter(rm_assessments=['Bad'], del_qc_var=False)\n",
    "\n",
    "ds_acsm.clean.cleanup()\n",
    "ds_acsm = act.qc.arm.add_dqr_to_qc(ds_acsm)\n",
    "ds_acsm.qcfilter.datafilter(rm_assessments=['Bad'], del_qc_var=False)\n",
    "\n",
    "ds_sp2.clean.cleanup()\n",
    "ds_sp2 = act.qc.arm.add_dqr_to_qc(ds_sp2)\n",
    "ds_sp2.qcfilter.datafilter(rm_assessments=['Bad'], del_qc_var=False)\n",
    "\n",
    "ds_psap.clean.cleanup()\n",
    "ds_psap = act.qc.arm.add_dqr_to_qc(ds_psap)\n",
    "ds_psap.qcfilter.datafilter(rm_assessments=['Bad'], del_qc_var=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af7d3e43-1607-4037-961b-8c09ea308f55",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# We can pass a dictionary to the display objects with multiple datasets\n",
    "# So let's plot all this up!\n",
    "display = act.plotting.TimeSeriesDisplay({'aps': ds_aps, 'mpl': ds_mpl, 'smps': ds_smps, 'acsm': ds_acsm, 'psap': ds_psap, 'sp2': ds_sp2, 'met': ds_met},\n",
    "                                         subplot_shape=(7,), figsize=(10,18))\n",
    "\n",
    "# MPL Plot\n",
    "# Variable names of interest linear_depol_ratio, linear_depol_snr, backscatter_snr\n",
    "display.plot('linear_depol_ratio', dsname='mpl', subplot_index=(0,), cvd_friendly=True)\n",
    "display.set_yrng([0, 3], subplot_index=(0,))\n",
    "\n",
    "# APS Plot\n",
    "display.plot('total_N_conc', dsname='aps', subplot_index=(1,))\n",
    "display.day_night_background(dsname='aps', subplot_index=(1,))\n",
    "\n",
    "# SMPS Plot\n",
    "display.plot('total_N_conc', dsname='smps', subplot_index=(2,))\n",
    "display.day_night_background(dsname='smps', subplot_index=(2,))\n",
    "\n",
    "# ACSM plot\n",
    "display.plot('sulfate', dsname='acsm', subplot_index=(3,), label='sulfate')\n",
    "display.plot('nitrate', dsname='acsm', subplot_index=(3,), label='nitrate')\n",
    "display.plot('ammonium', dsname='acsm', subplot_index=(3,), label='ammonium')\n",
    "display.plot('chloride', dsname='acsm', subplot_index=(3,), label='chloride')\n",
    "display.plot('total_organics', dsname='acsm', subplot_index=(3,), label='total_organics')\n",
    "\n",
    "display.day_night_background(dsname='acsm', subplot_index=(3,))\n",
    "\n",
    "# SP2 Plot\n",
    "display.plot('sp2_rbc_conc', dsname='sp2', subplot_index=(4,))\n",
    "display.day_night_background(dsname='sp2', subplot_index=(4,))\n",
    "\n",
    "# PSAP Plot\n",
    "display.plot('Ba_B_Weiss', dsname='psap', subplot_index=(5,))\n",
    "display.day_night_background(dsname='psap', subplot_index=(5,))\n",
    "\n",
    "# MET Wind Plot\n",
    "display.plot('wdir_vec_mean', dsname='met', subplot_index=(6,), linestyle='', marker='.')\n",
    "#display.plot('wspd_vec_mean', dsname='met', subplot_index=(6,), marker=None, secondary_y=True, color='orange', label='Wind Speed')\n",
    "display.day_night_background(dsname='met', subplot_index=(6,))\n",
    "\n",
    "plt.subplots_adjust(hspace=0.3)\n",
    "plt.legend()\n",
    "plt.savefig('./images/output.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1b979c8-f9c5-4420-8845-9b73b70fbf1b",
   "metadata": {},
   "source": [
    "#### Check this out\n",
    "[Floutsi et al 2023](https://amt.copernicus.org/articles/16/2353/2023/) has a great paper that shows how different aerosols can have different ratios.  In [Figure 3](https://amt.copernicus.org/articles/16/2353/2023/amt-16-2353-2023-f03-web.png), you can see that Sharan dust has a slightly different signal than Central Asian dust."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81e99004-3e0b-474e-b016-8e5e79a161d0",
   "metadata": {},
   "source": [
    "### Try more data!\n",
    "\n",
    "Go back and adjust the start/end date around this period and see how things change!  You can also try adding in some of the other instruments like the OPC (houaosopcM1.b1: calc_dust_weight or total_N_conc) or the NEPHDRY (houaosnephdry1mM1.b1: Bs_B_Dry_Neph3W).  You can find a full list of data products on the ARM TRACER [website](https://www.arm.gov/research/campaigns/amf2021tracer#:~:text=Order%20Data-,HOU%20DATA%20SOURCES,-NAME)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcef61c0-becf-4f59-a92a-7707997877ae",
   "metadata": {},
   "source": [
    "## Advanced Visualizations and Processing\n",
    "Let's try and dive into the data a little more and see if there are any patterns in the data based on the direction of the wind.\n",
    "\n",
    "### Data Roses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53326c0d-df46-4ce6-b0a1-a87e6bd90641",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# We already should have the data loaded up so let's explore with some data roses\n",
    "# First we need to combine data and to do that, we need to get it on the same time grid\n",
    "ds_combined = xr.merge([ds_met.resample(time='30min').nearest(), ds_acsm.resample(time='30min').nearest()], compat='override')\n",
    "\n",
    "# Plot out the data rose using the WindRose display object\n",
    "display = act.plotting.WindRoseDisplay(ds_combined)\n",
    "display.plot_data('wdir_vec_mean', 'wspd_vec_mean', 'sulfate', num_dirs=15, plot_type='line', line_plot_calc='mean')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fe0eb42-34d3-4e58-ac3b-49ac6c0453d8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# First we need to combine data and to do that, we need to get it on the same time grid\n",
    "ds_combined = xr.merge([ds_met.resample(time='1min').nearest(), ds_sp2.resample(time='1min').nearest()], compat='override')\n",
    "\n",
    "# Plot out the data rose using the WindRose display object\n",
    "display = act.plotting.WindRoseDisplay(ds_combined)\n",
    "# Let's try a different type of data rose that will show the mean Black Carbon Concentration\n",
    "# depending on wind direction and speed\n",
    "display.plot_data('wdir_vec_mean', 'wspd_vec_mean', 'sp2_rbc_conc', num_dirs=15, plot_type='contour', contour_type='mean')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa7ae058-c48d-48da-bbe9-fb2ce1c6c26f",
   "metadata": {},
   "source": [
    "### Checkout the area\n",
    "The AMF was deployed at [La Porte Municipal Airport](https://www.google.com/maps/place/Airport+Blvd/@29.6652378,-95.0466689,9165m/data=!3m1!1e3!4m7!3m6!1s0x863f6020e5e0ea21:0x792ee34f8eaac3e8!4b1!8m2!3d29.6663473!4d-95.0578571!16s%2Fg%2F1wbf_smp?entry=ttu).  Check out the google map and see if this mapes sense!\n",
    "\n",
    "### Now back to the dust event!\n",
    "Let's see if we can look at the lowest layers in the MPL data and see how they compare to the APS data.  For this, we can use simple matplotlib plots to investigate as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1560e9f0-a398-4a05-a438-9e0216ca52b1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Let's take a slice of the data to figure out where the signal\n",
    "# may be showing up in the profile\n",
    "data = ds_mpl['linear_depol_ratio'].values[:, 0:10]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(14,6))\n",
    "ax.plot(ds_mpl['time'], data)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dadae53c-a37d-424b-9bb2-5c2def73ec86",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Let's do it again and add in the APS data\n",
    "# Other variables to try: backscatter_snr, linear_depol_snr, linear_depol_ratio\n",
    "data = ds_mpl['linear_depol_ratio'].values[:, 4]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12,6))\n",
    "ax.plot(ds_mpl['time'], data, label='MPL LDR')\n",
    "\n",
    "ax2 = ax.twinx()\n",
    "ax2.plot(ds_aps['time'], ds_aps['total_N_conc'], color='purple', label='APS')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e2d172c-c383-4d4f-917a-59338a82d00d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Let's try a scatter plot using matplotlib\n",
    "# First we need the data on the same time scale\n",
    "ds_mpl2 = ds_mpl.resample(time='60s').nearest()\n",
    "ds_aps2 = ds_aps.resample(time='60s').nearest()\n",
    "\n",
    "ds = xr.merge([ds_mpl2, ds_aps2], compat='override')\n",
    "\n",
    "# Other variables to try: backscatter_snr, linear_depol_snr, linear_depol_ratio\n",
    "data = ds['linear_depol_snr'].values[:, 4]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 8))\n",
    "ax.scatter(data, ds['total_N_conc']) # , c=ds['linear_depol_snr'].values[:, 4])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4178c3dd-4a15-4e05-839c-b5fdf96c0efd",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec6021b0",
   "metadata": {},
   "source": [
    "### Excercises\n",
    "Here are a couple additional excercises that users can work through on their own time."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbc5bd72",
   "metadata": {},
   "source": [
    "#### Time series of ratios of fine particles vs. coarse particles using SMPS and APS data\n",
    "Using the SMPS (ds_smps) and APS (ds_aps) data, can you calculate a ratio of fine to coarse particles?  Remember, that this data would need to be on the same timestep."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "830f729b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, let's resample to 1-minute averages.  Note, that this may take some time with this much data.\n",
    "# The other option is to grab the nearest value\n",
    "ds_smps_1m = ds_smps.resample(time='1min').nearest()\n",
    "ds_aps_1m = ds_aps.resample(time='1min').nearest()\n",
    "\n",
    "# But wait, the concentration variables have the same name!  We need to rename\n",
    "ds_aps_1m = ds_aps_1m.rename_vars({'total_N_conc': 'total_N_conc_aps'})\n",
    "ds_merge = xr.merge([ds_smps_1m, ds_aps_1m], compat='override')\n",
    "\n",
    "# Next, let's calculate the ratio\n",
    "ratio = ds_merge['total_N_conc'].values / ds_merge['total_N_conc_aps'].values\n",
    "\n",
    "# Then we can add it back to the Dataset as a variable\n",
    "atts = {'units': '1', 'long_name': 'Ration of SMPS to APS total_N_conc'}\n",
    "da = xr.DataArray(ratio, coords=ds_merge['total_N_conc'].coords, dims=ds_merge['total_N_conc'].dims, attrs=atts)\n",
    "ds_merge['smps_aps_ratio'] = da\n",
    "\n",
    "# And then plot it wil ACT!\n",
    "display = act.plotting.TimeSeriesDisplay(ds_merge, figsize=(12, 10), subplot_shape=(3,))\n",
    "display.plot('smps_aps_ratio', subplot_index=(0,))\n",
    "display.plot('total_N_conc', subplot_index=(1,))\n",
    "display.plot('total_N_conc_aps', subplot_index=(2,))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dea39c3",
   "metadata": {},
   "source": [
    "#### Compare Chemical Composition Between Events\n",
    "\n",
    "Create a pie chart of aerosol composition on July 12, 13, and compare it to that from 17-18.  This will require subsetting the ds_acsm data.\n",
    "\n",
    "Studies have shown that Dust events can contribute to sulfate formation through the interaction of mineral dust particles with sulfur dioxide (SO2) and other atmospheric gases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70f8a3c7-8394-4fd7-8d93-c255b113f7b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# First subset the data for both sets\n",
    "ds_acsm_12 = ds_acsm.sel(time=slice('2022-07-12', '2022-07-13'))\n",
    "ds_acsm_17 = ds_acsm.sel(time=slice('2022-07-17', '2022-07-18'))\n",
    "\n",
    "# Run through and get data\n",
    "fields = ['total_organics', 'sulfate', 'nitrate', 'ammonium', 'chloride']\n",
    "\n",
    "# We also want to provide some keyword arguments to avoid invalid data such\n",
    "# as negative values.\n",
    "threshold = 0.0\n",
    "fill_value = 0.0\n",
    "\n",
    "# Create a DistributionDisplay object to compare fields\n",
    "display = act.plotting.DistributionDisplay({'acsm_12': ds_acsm_12, 'acsm_17': ds_acsm_17}, \n",
    "                                           subplot_shape=(1,2), figsize=(10,14))\n",
    "\n",
    "# We can set one of the slices to explode and give it a nice shadow.\n",
    "explode = (0, 0.1, 0, 0, 0)\n",
    "shadow = True\n",
    "\n",
    "# Create a pie chart using the fields list. The percentages of the\n",
    "# fields will be calculated using act.utils.calculate_percentages.\n",
    "display.plot_pie_chart(\n",
    "    fields,\n",
    "    dsname='acsm_12',\n",
    "    threshold=threshold,\n",
    "    fill_value=fill_value,\n",
    "    explode=explode,\n",
    "    shadow=True,\n",
    "    subplot_index=(0,0)\n",
    ")\n",
    "display.plot_pie_chart(\n",
    "    fields,\n",
    "    dsname='acsm_17',\n",
    "    threshold=threshold,\n",
    "    fill_value=fill_value,\n",
    "    explode=explode,\n",
    "    shadow=True,\n",
    "    subplot_index=(0,1)\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc0c5efb-0a18-4fdc-a697-1d299b77b0e0",
   "metadata": {},
   "source": [
    "#### Chemical Composition Distribution\n",
    "\n",
    "The PI charts can provide some information on the distribution but we can get more information by looking at the distribution using violin plots.  Does the dust event still show high sulfates in this view?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee7243eb-b2ba-48d0-aac4-8c48654a6c71",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a DistributionDisplay object to compare fields\n",
    "display = act.plotting.DistributionDisplay({'acsm_12': ds_acsm_12, 'acsm_17': ds_acsm_17},\n",
    "                                           subplot_shape=(1,2), figsize=(12,6))\n",
    "\n",
    "# Compare the distributions of each composition type\n",
    "display.plot_violin('total_organics', positions=[1.0], dsname='acsm_12', subplot_index=(0, 0))\n",
    "display.plot_violin('sulfate', positions=[2.0], dsname='acsm_12', subplot_index=(0, 0))\n",
    "display.plot_violin('ammonium', positions=[3.0], dsname='acsm_12', subplot_index=(0, 0))\n",
    "display.plot_violin('nitrate', positions=[4.0], dsname='acsm_12', subplot_index=(0, 0))\n",
    "display.plot_violin('chloride', positions=[5.0], dsname='acsm_12', subplot_index=(0, 0))\n",
    "\n",
    "display.plot_violin('total_organics', positions=[1.0], dsname='acsm_17', subplot_index=(0, 1))\n",
    "display.plot_violin('sulfate', positions=[2.0], dsname='acsm_17', subplot_index=(0, 1))\n",
    "display.plot_violin('ammonium', positions=[3.0], dsname='acsm_17', subplot_index=(0, 1))\n",
    "display.plot_violin('nitrate', positions=[4.0], dsname='acsm_17', subplot_index=(0, 1))\n",
    "display.plot_violin('chloride', positions=[5.0], dsname='acsm_17', subplot_index=(0, 1))\n",
    "\n",
    "# Update the tick information\n",
    "ticks = ['', 'Total Organics', 'Sulfate', 'Ammonium', 'Nitrate', 'Chloride', '']\n",
    "display.axes[0, 0].set_xticks([0.5, 1, 2, 3, 4, 5, 5.5])\n",
    "display.axes[0, 1].set_xticks([0.5, 1, 2, 3, 4, 5, 5.5])\n",
    "display.axes[0, 0].set_xticklabels(ticks)\n",
    "display.axes[0, 1].set_xticklabels(ticks)\n",
    "display.axes[0, 0].set_ylim([0, 3.5])\n",
    "display.axes[0, 1].set_ylim([0, 3.5])\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac2ef64d-14ba-4bee-a578-46137c7ebb03",
   "metadata": {},
   "source": [
    "## Doppler Lidar\n",
    "Out of interest, let us take a look at the doppler lidar data.\n",
    "\n",
    "We can used both the fixed pointing mode (FPT) and plan position indicator (PPI) mode datasets to plot signal to noise ration\n",
    "with an overlay of wind barbs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec38347a-1973-41aa-a0ee-aa08ee4dd7bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download the data as before\n",
    "# Read an ARM TRACER Doppler lidar FPT dataset\n",
    "datastream = 'houdlfptM1.b1'\n",
    "startdate = '2022-07-16T00:00:00'\n",
    "enddate = '2022-07-16T23:50:00'\n",
    "result_dl_fpt = glob.glob('./' + datastream +'/*')\n",
    "if len(result_dl_fpt) == 0:\n",
    "    result_dl_fpt = act.discovery.download_arm_data(username, token, datastream, startdate, enddate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7127d4d-a7cc-4463-9e9e-3f3dd351788c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets read all the files into one timeseries dataset.\n",
    "ds_dlfpt = act.io.read_arm_netcdf(result_dl_fpt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe9736f0-81ec-4bd1-aa57-bcdc64a18e96",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change intensity to signal to noise ratio\n",
    "data = ds_dlfpt['intensity'].data\n",
    "ds_dlfpt['intensity'].data = 1. - np.log10(data)\n",
    "ds_dlfpt['intensity'].attrs['units'] = '1 - LOG10 ' + ds_dlfpt['intensity'].attrs['long_name']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60e79a81-6392-47a3-8af1-126caf9f4048",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Next, let's read in the ARM TRACER Doppler lidar PPI dataset\n",
    "datastream = 'houdlppiM1.b1'\n",
    "result_dl_ppi = glob.glob('./' + datastream +'/*')\n",
    "if len(result_dl_ppi) == 0:\n",
    "    result_dl_ppi = act.discovery.download_arm_data(username, token, datastream, startdate, enddate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d20c8964-4cf8-4d58-ae4a-7f37b1d70262",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instead of reading in all files to a dataset, lets calculate the winds for each file and\n",
    "# merge them into one dataset to use for plotting the wind barbs.\n",
    "multi_ds = []\n",
    "for file in sorted(result_dl_ppi):\n",
    "    ds = act.io.arm.read_arm_netcdf(file)\n",
    "    # Calculate the winds for each houdlppi dataset.\n",
    "    wind_ds = act.retrievals.compute_winds_from_ppi(\n",
    "        ds, remove_all_missing=False, snr_threshold=0.008\n",
    "    )\n",
    "    multi_ds.append(wind_ds)\n",
    "\n",
    "wind_ds = xr.merge(multi_ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d1cebb5-c388-456f-9719-f1c0889cfc20",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's create a display object containing both the houdlfpt and calculated winds dataset\n",
    "display = act.plotting.TimeSeriesDisplay(\n",
    "    {\n",
    "        \"HOU DLPPI Computed Winds Over DL FPT Backscatter\": wind_ds,\n",
    "        \"houdlfptM1.b1\": ds_dlfpt,\n",
    "    },\n",
    "    figsize=(20, 10),\n",
    ")\n",
    "\n",
    "# Plot the wind barbs overlayed on the signal to noise ratio\n",
    "display.plot(\n",
    "    'intensity', dsname='houdlfptM1.b1', cmap='NWSRef', vmax=1.0, vmin=0.0\n",
    ")\n",
    "display.plot_barbs_from_spd_dir(\n",
    "    'wind_speed', 'wind_direction', dsname='HOU DLPPI Computed Winds Over DL FPT Backscatter', invert_y_axis=False, barbcolor='white', barb_step_y=2\n",
    ")\n",
    "\n",
    "# Update the y-limits to show plotted winds\n",
    "display.axes[0].set_ylim([0, 2500])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7de4682-fb16-48ac-bc5a-29a8385dd29c",
   "metadata": {},
   "source": [
    "### What's Happening Here?\n",
    "Take a look at the [DQ Office Plots](https://plot.adc.arm.gov/PLOTS/hou/houmet/20220722/houmetM1.b1.meteogram.20220722.png) of surface meteorological variables for the day. In particular, check out the wind direction and see if that matches up with what you're seeing in the lowest level of the plots."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd37d2fc",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1160da20-2f54-43c0-aeea-d7a306d37f6d",
   "metadata": {},
   "source": [
    "## Additional Questions for the User to Explore\n",
    "1. Can you improve these matplotlib plots to make them more readable?\n",
    "1. Are there other variables that you would be interested in for the data roses?\n",
    "1. We looked at relationships between the MPL and APS, are there other sets of instruments that might have relationships? (SP2 and PSAP?)\n",
    "\n",
    "#### More Advanced\n",
    "1. Can we use a machine learning toolkit like scikit-learn to cluster the MPL/APS data?\n",
    "1. Could we get relationships between those clusters or a classification of what it is (smoke, dust, pollution)?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a810279c-09ba-49e7-b517-e791c350883b",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0eaee516-805e-46ed-94e9-1d115bcebb6c",
   "metadata": {},
   "source": [
    "## Next Steps\n",
    "This notebook showed users how they can explore multiple different datasets using ACT.  The next steps are for the users to continue exploring and diving into the data and the questions posed at the end!\n",
    "\n",
    "### Data Used in this Notebook\n",
    "Ermold, B., & Flynn, C. Particle Soot Absorption Photometer (AOSPSAP3W1M). Atmospheric Radiation Measurement (ARM) User Facility. https://doi.org/10.5439/1225037\n",
    "\n",
    "Jackson, R., & Sedlacek, A. Single Particle Soot Photometer (AOSSP2BC60S). Atmospheric Radiation Measurement (ARM) User Facility. https://doi.org/10.5439/1807910\n",
    "\n",
    "Kuang, C., & Singh, A. Aerodynamic Particle Sizer (AOSAPS). Atmospheric Radiation Measurement (ARM) User Facility. https://doi.org/10.5439/1407135\n",
    "\n",
    "Kyrouac, J., & Shi, Y. Surface Meteorological Instrumentation (MET). Atmospheric Radiation Measurement (ARM) User Facility. https://doi.org/10.5439/1786358\n",
    "\n",
    "Muradyan, P., Cromwell, E., Koontz, A., & Coulter, R. Micropulse Lidar (MPLPOLFS). Atmospheric Radiation Measurement (ARM) User Facility. https://doi.org/10.5439/1320657\n",
    "\n",
    "Newsom, R., Shi, Y., & Krishnamurthy, R. Doppler Lidar (DLFPT). Atmospheric Radiation Measurement (ARM) User Facility. https://doi.org/10.5439/1025185\n",
    "\n",
    "Newsom, R., Shi, Y., & Krishnamurthy, R. Doppler Lidar (DLPPI). Atmospheric Radiation Measurement (ARM) User Facility. https://doi.org/10.5439/1025186\n",
    "\n",
    "Sivaraman, C., Flynn, D., Riihimaki, L., & Comstock, J. Cloud mask from Micropulse Lidar (30SMPLCMASK1ZWANG). Atmospheric Radiation Measurement (ARM)\n",
    "User Facility. https://doi.org/10.5439/1508389\n",
    "\n",
    "Zawadowicz, M., & Howie, J. Aerosol Chemical Speciation Monitor (AOSACSM). Atmospheric Radiation Measurement (ARM) User Facility. https://doi.org/10.5439/1762267"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}