01_Masking.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "toc_visible": true,
      "authorship_tag": "ABX9TyNgeV1x3cjoszIpFZFGOYNJ",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/aamatthews/2022UWCourse/blob/main/01_Masking.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Setup"
      ],
      "metadata": {
        "id": "j88FkfBVi2KS"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "At the start of each notebook, we will need to grab the data again, so here are those lines!"
      ],
      "metadata": {
        "id": "QppnvRs7NXH3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!wget https://engineering.arm.gov/~amatthews/hou_data_sample/20220701.zip"
      ],
      "metadata": {
        "id": "iz_6mTVkNiui"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!unzip 20220701.zip"
      ],
      "metadata": {
        "id": "otoHEBMwNklD"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import netCDF4 as nc\n",
        "import xarray as xr\n",
        "import pandas as pd\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline"
      ],
      "metadata": {
        "id": "pEZWwmQMPME8"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Introduction"
      ],
      "metadata": {
        "id": "oLwesHOHjWmS"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "As you may have noticed in the data we plotted previously, there was a lot of data that appeared to not be weather. \n",
        "\n",
        "In this next part, you will learn methods to hide or 'mask out' the non-meteorological data."
      ],
      "metadata": {
        "id": "nVCYlQFbOJRG"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "First, lets read in a data file, get a list of the variables, and plot reflectivity."
      ],
      "metadata": {
        "id": "WCCslUI8jLM9"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "ppi = xr.open_dataset('/content/20220701/houcsapr2cfrS2.a1.20220701.192603.nc', decode_times = False)\n",
        "ppi['time'] = pd.to_datetime(ppi['base_time'].values + ppi['time_offset'].values, unit='s')"
      ],
      "metadata": {
        "id": "xAy8RrUwNooy"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "list(ppi.variables.keys())"
      ],
      "metadata": {
        "id": "4CwTtJyWQTED"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "sweep_start = list(map(int,(ppi['sweep_start_ray_index'].data)))\n",
        "sweep_end = list(map(int,(ppi['sweep_end_ray_index'].data)))\n",
        "azi = ppi['azimuth'].data\n",
        "elv = ppi['elevation'].data\n",
        "r = ppi['range'].data\n",
        "\n",
        "y = np.outer(r, np.cos((azi)*np.pi/180.0))\n",
        "x = np.outer(r, np.sin((azi)*np.pi/180.0))\n",
        "sweep_num = 0\n",
        "\n",
        "\n",
        "sweep_num = 1\n",
        "#     print(sweep_num)\n",
        "plt.figure(figsize=(10,7.5))\n",
        "\n",
        "plt.pcolormesh(x[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                y[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                ppi['reflectivity'].data[sweep_start[sweep_num]:sweep_end[sweep_num],:], \n",
        "                vmin=-60, vmax=40, cmap='jet')\n",
        "clb = plt.colorbar()\n",
        "clb.set_label('Zh')\n",
        "plt.xlabel('E-W Distance from Radar (m)')\n",
        "plt.ylabel('N-S Distance from Radar (m)')\n",
        "elv_start = elv[sweep_start[sweep_num]+50]\n",
        "plt.title('Elevation: %f' % elv_start)"
      ],
      "metadata": {
        "id": "Zgqh_ymRPKdn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Masking"
      ],
      "metadata": {
        "id": "td427QYdjgAw"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "One variable we often use to find the meteorological signal is Signal to Noise Ratio, as when this is higher there is less noise in the return."
      ],
      "metadata": {
        "id": "14xB1RgwRJRy"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plt.figure(figsize=(10,7.5))\n",
        "\n",
        "plt.pcolormesh(x[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                y[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                ppi['signal_to_noise_ratio_copolar_h'].data[sweep_start[sweep_num]:sweep_end[sweep_num],:], \n",
        "                vmin=-60, vmax=40, cmap='jet')\n",
        "clb = plt.colorbar()\n",
        "clb.set_label('SNR_h')\n",
        "plt.xlabel('E-W Distance from Radar (m)')\n",
        "plt.ylabel('N-S Distance from Radar (m)')\n",
        "elv_start = elv[sweep_start[sweep_num]+50]\n",
        "plt.title('Elevation: %f' % elv_start)"
      ],
      "metadata": {
        "id": "fTU-9bP6QDPi"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's also take a look at the noise floor using a histogram."
      ],
      "metadata": {
        "id": "ENXrP_6nkSbR"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plt.hist(ppi['signal_to_noise_ratio_copolar_h'].data.flatten(), bins=100);"
      ],
      "metadata": {
        "id": "88J7dw4_jy2Z"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Based on these plots, what do you think would be a good SNR value to start with as our threshold?"
      ],
      "metadata": {
        "id": "Q9iKSm37RV8I"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "snr_thres = -5"
      ],
      "metadata": {
        "id": "7iLIrWwMRFv7"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "mask_snr = np.zeros_like(ppi['signal_to_noise_ratio_copolar_h'].data)"
      ],
      "metadata": {
        "id": "REpiJz-HQR6a"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "mask_snr[ppi['signal_to_noise_ratio_copolar_h'].data > snr_thres] = True"
      ],
      "metadata": {
        "id": "gO9A1_kRRCf0"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "mask_snr"
      ],
      "metadata": {
        "id": "Ge8_zLhQRugA"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plt.figure(figsize=(10,7.5))\n",
        "\n",
        "plt.pcolormesh(x[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                y[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                mask_snr[sweep_start[sweep_num]:sweep_end[sweep_num],:], \n",
        "                vmin=0, vmax=1)\n",
        "clb = plt.colorbar()\n",
        "clb.set_label('SNR_mask')\n",
        "plt.xlabel('E-W Distance from Radar (m)')\n",
        "plt.ylabel('N-S Distance from Radar (m)')\n",
        "elv_start = elv[sweep_start[sweep_num]+50]\n",
        "plt.title('Elevation: %f' % elv_start)"
      ],
      "metadata": {
        "id": "OFkvJuzfRDbr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "As you can see, this only removes the true noise values from the data, but not all the returns that are non-meteorological. Let's go back and edit the above threshold to be higher. \n",
        "\n",
        "One issue with an SNR mask, and many masking methods in general, is that you have to balance removing background with not removing too much cloud. If you make the SNR threshold too high, then you can cut off the edges of the cloud. It all depends on what you are interested in for your research."
      ],
      "metadata": {
        "id": "yIgPoPu1kbz3"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's plot RhoHV next and see if it can help us filter the data. \n"
      ],
      "metadata": {
        "id": "WUcxuKZJR9HJ"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "When working with new variables, it can be helpful to print the minumum and maximum values so you know what to set as your plot mins and maxes."
      ],
      "metadata": {
        "id": "h708jh2Vk5My"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "np.nanmin(ppi['copol_correlation_coeff'])"
      ],
      "metadata": {
        "id": "q8UKRp4kSXIt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "np.nanmax(ppi['copol_correlation_coeff'])"
      ],
      "metadata": {
        "id": "MkyVpeZRSTou"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plt.figure(figsize=(10,7.5))\n",
        "\n",
        "plt.pcolormesh(x[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                y[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                ppi['copol_correlation_coeff'].data[sweep_start[sweep_num]:sweep_end[sweep_num],:], \n",
        "                vmin=0, vmax=1, cmap='jet')\n",
        "clb = plt.colorbar()\n",
        "clb.set_label('RhoHV')\n",
        "plt.xlabel('E-W Distance from Radar (m)')\n",
        "plt.ylabel('N-S Distance from Radar (m)')\n",
        "elv_start = elv[sweep_start[sweep_num]+50]\n",
        "plt.title('Elevation: %f' % elv_start)"
      ],
      "metadata": {
        "id": "ComhAcguR3ud"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "A RhoHV value for a metorological signal should be close to 1, as it is the correlation between the H and V returns. A rougher surface (like the ground, or a bird) will have a lower RhoHV value."
      ],
      "metadata": {
        "id": "ZcLXKuhclLXf"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's filter with a threshold of 0.95 to see how that looks, and we can increase that if we would like."
      ],
      "metadata": {
        "id": "zAgJ-oHVSgl5"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "rhohv_thres = 0.95"
      ],
      "metadata": {
        "id": "OIVz3sjsSPLn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "mask_rhohv = np.zeros_like(ppi['copol_correlation_coeff'])"
      ],
      "metadata": {
        "id": "qQxqHNyySvLK"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "mask_rhohv[ppi['copol_correlation_coeff'].data > rhohv_thres] = True"
      ],
      "metadata": {
        "id": "2J4cPN37Sy2J"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plt.figure(figsize=(10,7.5))\n",
        "\n",
        "plt.pcolormesh(x[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                y[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                mask_rhohv[sweep_start[sweep_num]:sweep_end[sweep_num],:], \n",
        "                vmin=0, vmax=1)\n",
        "clb = plt.colorbar()\n",
        "clb.set_label('SNR_mask')\n",
        "plt.xlabel('E-W Distance from Radar (m)')\n",
        "plt.ylabel('N-S Distance from Radar (m)')\n",
        "elv_start = elv[sweep_start[sweep_num]+50]\n",
        "plt.title('Elevation: %f' % elv_start)"
      ],
      "metadata": {
        "id": "m_-YumYwSzyZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Wow! This filters out much more than the SNR mask did! What happens when we combine them?"
      ],
      "metadata": {
        "id": "nWzRMAMNTF4I"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "mask_all = np.logical_and(mask_snr, mask_rhohv)"
      ],
      "metadata": {
        "id": "iMZKCfX7TES6"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plt.figure(figsize=[15,12.5])\n",
        "\n",
        "plt.subplot(221)\n",
        "\n",
        "plt.pcolormesh(x[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                y[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                mask_snr[sweep_start[sweep_num]:sweep_end[sweep_num],:], \n",
        "                vmin=0, vmax=1)\n",
        "clb = plt.colorbar()\n",
        "clb.set_label('Mask SNR')\n",
        "plt.xlabel('E-W Distance from Radar (m)')\n",
        "plt.ylabel('N-S Distance from Radar (m)')\n",
        "elv_start = elv[sweep_start[sweep_num]+50]\n",
        "plt.title('SNR Mask \\n Elevation: %f' % elv_start)\n",
        "\n",
        "plt.subplot(222)\n",
        "plt.pcolormesh(x[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                y[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                mask_rhohv[sweep_start[sweep_num]:sweep_end[sweep_num],:], \n",
        "                vmin=0, vmax=1)\n",
        "clb = plt.colorbar()\n",
        "clb.set_label('Mask Rhohv')\n",
        "plt.xlabel('E-W Distance from Radar (m)')\n",
        "plt.ylabel('N-S Distance from Radar (m)')\n",
        "elv_start = elv[sweep_start[sweep_num]+50]\n",
        "plt.title('RhoHV Mask \\n Elevation: %f' % elv_start)\n",
        "\n",
        "plt.subplot(223)\n",
        "plt.pcolormesh(x[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                y[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                mask_all[sweep_start[sweep_num]:sweep_end[sweep_num],:], \n",
        "                vmin=0, vmax=1)\n",
        "clb = plt.colorbar()\n",
        "clb.set_label('Mask')\n",
        "plt.xlabel('E-W Distance from Radar (m)')\n",
        "plt.ylabel('N-S Distance from Radar (m)')\n",
        "elv_start = elv[sweep_start[sweep_num]+50]\n",
        "plt.title('Combined Mask \\n Elevation: %f' % elv_start);"
      ],
      "metadata": {
        "id": "sjUUHAh0TR9b"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "So, what does this look like when applied to the reflectivity?"
      ],
      "metadata": {
        "id": "xc2K5t1bgDWP"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "z_masked = np.ma.masked_where(mask_all == False, ppi['reflectivity'])"
      ],
      "metadata": {
        "id": "241S7mjUTVFr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plt.figure(figsize=(14,5))\n",
        "\n",
        "plt.subplot(121)\n",
        "plt.pcolormesh(x[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                y[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                ppi['reflectivity'].data[sweep_start[sweep_num]:sweep_end[sweep_num],:], \n",
        "                vmin=-60, vmax=40, cmap='jet')\n",
        "clb = plt.colorbar()\n",
        "clb.set_label('Zh')\n",
        "plt.xlabel('E-W Distance from Radar (m)')\n",
        "plt.ylabel('N-S Distance from Radar (m)')\n",
        "elv_start = elv[sweep_start[sweep_num]+50]\n",
        "plt.title('Raw Reflectivity \\n Elevation: %f' % elv_start)\n",
        "\n",
        "\n",
        "plt.subplot(122)\n",
        "plt.pcolormesh(x[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                y[:,sweep_start[sweep_num]:sweep_end[sweep_num]].T,\n",
        "                z_masked[sweep_start[sweep_num]:sweep_end[sweep_num],:], \n",
        "                vmin=-60, vmax=40, cmap='jet')\n",
        "clb = plt.colorbar()\n",
        "clb.set_label('Zh')\n",
        "plt.xlabel('E-W Distance from Radar (m)')\n",
        "plt.ylabel('N-S Distance from Radar (m)')\n",
        "elv_start = elv[sweep_start[sweep_num]+50]\n",
        "plt.title('Masked Reflectivity \\n Elevation: %f' % elv_start);"
      ],
      "metadata": {
        "id": "qO8Ty5tBgUcR"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "As you can see, this removes a good portion of the non-meteorological echo, but not all! There are still some extra features present, both likely clutter signals and possible second trip echos (an echo from a given pulse that is not received until after the transmission of the next pulse) \n",
        "\n"
      ],
      "metadata": {
        "id": "r7PMsrVpg9kD"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Exercise\n",
        "1. Are there any additional variables you think could be included in the mask? Test them out!\n",
        "\n",
        "2. Will the full mask you came up with work for the X and Ka bands as well? Or will the numbers we used here need to be modified? Test this on them!"
      ],
      "metadata": {
        "id": "_wZIE9rcmv61"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "CgqUtffwgcF0"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}