tutorials/opac7irowg9_workshop.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to AWS cloud computing and GNSS RO Data in the AWS Open Data Registry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Authors:** Stephen Leroy and Amy McVey\n",
    "\n",
    "**Date:** September 13, 2022\n",
    "\n",
    "This notebook is composed for instruction and experimentation at the workshop to \n",
    "introduce cloud computing and how it can be applied to GNSS radio occultation \n",
    "data in the AWS Open Data Registry. Feel free to write any of your own code \n",
    "based on this notebook in your EC2 instance. Python is available in it for your \n",
    "convenience. \n",
    "\n",
    "- AWS portal to the RO repository in the Open Data Registry: \n",
    "    https://registry.opendata.aws/gnss-ro-opendata/\n",
    "\n",
    "- GitHub support material: \n",
    "    http://github.com/gnss-ro/aws-opendata/\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Portal to Open Data Registry of RO data with s3fs "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable \"s3\" is an s3fs object, which makes the contents of S3 buckets appear as if they were in a Linux-like file system. The AWS Open Data Registry is free for browse and download, so you need to inform AWS that no authentication (\"signature\") is necessary for access. That's why the signature_version is set to UNSIGNED. Also, the RO data is hosted in AWS region \"us-east-1\", and that needs to be declared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import s3fs\n",
    "from botocore import UNSIGNED\n",
    "\n",
    "AWS_region = \"us-east-1\"\n",
    "s3 = s3fs.S3FileSystem( \n",
    "        client_kwargs = { 'region_name': AWS_region },  \n",
    "        config_kwargs = { 'signature_version': UNSIGNED }\n",
    "        )   \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Attempt a few s3fs methods. First, list all RO missions contributed by UCAR. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "ucar_missions = sorted( [ os.path.split(p)[1] for p in s3.ls(\"gnss-ro-data/contributed/v1.1/ucar/\") ] )\n",
    "print( ucar_missions )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now list all files containing bending angle profiles\n",
    "as retrieved by UCAR for Metop satellites on January 1, 2020..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = sorted( s3.ls( \"gnss-ro-data/contributed/v1.1/ucar/metop/refractivityRetrieval/2020/01/01/\" ) )\n",
    "print( \"\\n\".join( files ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using a DynamoDB database"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The AWS Open Data Registry only allows accessible S3 storage of \n",
    "data; it does not allow publication of a database service. At home, \n",
    "you will need to use our GitHub utilities to construct your own \n",
    "private database based on the contents of \n",
    "s3://gnss-ro-data/dynamo/v1.1/export_subsets,\n",
    "which will take considerable time. At this workshop, we have made \n",
    "AER's own internal DynamoDB database of RO data available so you \n",
    "can learn how to query a DynamoDB database. Authentication has been \n",
    "enabled behind the scenes. \n",
    "\n",
    "First, the portal into the AWS Universe.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import boto3\n",
    "session = boto3.Session( region_name=AWS_region )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At home, it is possible that this won't work because the requisite security tokens giving you permission to access your own database are not in the computing environment. If the tokens are not available in the environment, then you will have to specify the current saml2aws profile you use for access to the database through the keyword *profile_name* in the call to boto3.Session. \n",
    "\n",
    "Commission the AWS service \"dynamodb\" by creating a \"dynamodb\" \n",
    "Python resource. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dynamodb_resource = session.resource( \"dynamodb\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Directly access a specific DynamoDB database, or \"table\". "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = dynamodb_resource.Table( \"gnss-ro-data-stagingv1_1\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every query requires a specific partition key, labeled as \n",
    "\"leo-ttt\" (receiver - dash - transmitter). This is how you \n",
    "define a \"key\" for the Metop-A receiver and occultations of \n",
    "GPS PRN 1 (\"G01\"): \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from boto3.dynamodb.conditions import Key\n",
    "partitionkey = Key('leo-ttt').eq( \"metopa-G01\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every query requires a sort key as well, but the sort key \n",
    "can be loosely defined, such as greater than or less than a value \n",
    "or between two values. Keep in mind that every combination of \n",
    "unique partition key and unique sort key must point to just one occultation \n",
    "sounding. Because an occultation is uniquely defined by the \n",
    "receiving satellite, the transmitting satellite, and the time \n",
    "of the sounding (precise to within a few minutes), all of those \n",
    "values must be contained in the partition key and sort key when \n",
    "combined. In our case, the partition key contains the receiver \n",
    "and the transmitter; therefore, the sort key must contain the time\n",
    "of the occultation sounding, and it does so under the label \"date-time\": "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sortkey = Key('date-time').between( \"2020-01-01-00-00\", \"2020-01-31-23-59\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Query the database. Notice how a compound key is formed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ret = table.query( KeyConditionExpression = partitionkey & sortkey )\n",
    "print( 'Query results found:', ret['Count'] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A filter on the results can be applied behind the scenes once \n",
    "AWS performs the query. That is done by filtering according to \n",
    "information in each RO item in the database. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from decimal import Decimal\n",
    "from boto3.dynamodb.conditions import Attr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter by region. Notice how compound filters are formed, using the \n",
    "binary \"and\" (\"&\") operator. The binary \"or\" operator (\"|\") and the \n",
    "logical \"not\" unary operator (\"~\") are also available. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filters = Attr('longitude').between( Decimal(120.0), Decimal(150.0) )\n",
    "filters = filters & Attr('latitude').between( Decimal(-35.0), Decimal(-15.0) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter for setting occultations only. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filters = filters & Attr(\"setting\").eq(\"True\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter for the availability of bending angle data as produced by \n",
    "the processing_center. Currently, *processing_center* can take on \n",
    "the values \"ucar\" or \"romsaf\". "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "processing_center = \"ucar\"\n",
    "filters = filters & Attr(f\"{processing_center}_refractivityRetrieval\").exists()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Query the database. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ret = table.query( KeyConditionExpression = partitionkey & sortkey,\n",
    "        FilterExpression = filters )\n",
    "print( 'Query results found:', ret['Count'] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The returned dictionary *ret* contains two keys of interest: \n",
    "*Count* and *Items*. *Count* is an integer count of the number \n",
    "of database items (occultation soundings) found by the query, and *Items* \n",
    "is a list of dictionaries, each dictionary containing all of the \n",
    "keys and information for one occultation sounding. \n",
    "\n",
    "Notice that these are occultations involving Metop-A and GPS 01 \n",
    "only. Now search for all Metop occultations in the same time period. \n",
    "A nested loop is necessary in order to specify all possible \n",
    "partition keys (over receiver and transmitter). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "receivers = [ \"metopa\", \"metopb\", \"metopc\" ]\n",
    "transmitters = [ f\"G{i:02d}\" for i in range(1,33) ]\n",
    "\n",
    "allitems = []       # Initialize the accumulative output list. \n",
    "\n",
    "for receiver in receivers:\n",
    "    for transmitter in transmitters:\n",
    "        partitionkey = Key('leo-ttt').eq( f\"{receiver}-{transmitter}\" )\n",
    "        ret = table.query(\n",
    "                KeyConditionExpression = partitionkey & sortkey,\n",
    "                FilterExpression = filters\n",
    "                )\n",
    "        allitems += ret['Items']    # Append to output list. \n",
    "\n",
    "print( f\"There were {len(allitems):d} soundings found.\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Congratulations! You have obtained a listing of all Metop occultation \n",
    "soundings, setting only, over Australia for the month of January, 2020. \n",
    "Now let's download all of those retrievals. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "local_paths = []\n",
    "os.makedirs('ro_data', exist_ok=True)\n",
    "\n",
    "for item in allitems: \n",
    "    refractivityRetrieval_file = item[f\"{processing_center}_refractivityRetrieval\"]\n",
    "    s3_path = f\"gnss-ro-data/{refractivityRetrieval_file}\"\n",
    "    local_path = os.path.join('ro_data',os.path.split( s3_path )[1])\n",
    "    local_paths.append( local_path )\n",
    "    if os.path.exists(local_path): continue\n",
    "    print( f\"Downloading {local_path}\" )\n",
    "    ret = s3.download( s3_path, local_path )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple analysis example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot all of the \"unoptimized\", ionospheric-corrected \n",
    "bending angles. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netCDF4 import Dataset\n",
    "import matplotlib.pyplot as plt \n",
    "from matplotlib.ticker import MultipleLocator\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the plot look nice. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axeslinewidth = 0.5 \n",
    "plt.rcParams.update( {\n",
    "  'font.family': \"Times New Roman\", \n",
    "  'font.size': 10,  \n",
    "  'font.weight': \"normal\",\n",
    "  'xtick.major.width': axeslinewidth,\n",
    "  'xtick.minor.width': axeslinewidth,\n",
    "  'ytick.major.width': axeslinewidth,\n",
    "  'ytick.minor.width': axeslinewidth,\n",
    "  'axes.linewidth': axeslinewidth } )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analysis program follows..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Define the x axis. \n",
    "\n",
    "fig = plt.figure( figsize=(6,6) )\n",
    "ax = fig.add_axes( [ 0.12, 0.12, 0.84, 0.86 ] )\n",
    "\n",
    "ax.set( \n",
    "    xlabel = \"Calibrated bending angle [mrad]\", \n",
    "    xticks = np.arange( 0.0, 40.01, 5.0 ), \n",
    "    xlim = ( 0.0, 40.0 ) \n",
    "    )\n",
    "ax.xaxis.set_minor_locator( MultipleLocator( 1.0 ) )\n",
    "\n",
    "#  Define the y axis. \n",
    "\n",
    "ax.set( \n",
    "    ylabel = \"Impact height [km]\", \n",
    "    yticks = np.arange( 0.0, 40.001, 10.0 ), \n",
    "    ylim = ( 0.0, 40.0 )\n",
    "    )\n",
    "ax.yaxis.set_minor_locator( MultipleLocator( 2.0 ) )\n",
    "\n",
    "#  Loop over files, converting impact parameter to impact height \n",
    "#  in km and bending angle to milliradians. \n",
    "\n",
    "for local_path in local_paths:\n",
    "\n",
    "    #  Read data. \n",
    "\n",
    "    d = Dataset( local_path, 'r' )\n",
    "    bending_angle = d.variables['bendingAngle'][:] * 1000.0\n",
    "    radius_of_curvature = d.variables['radiusOfCurvature'].getValue()\n",
    "    impact_height = ( d.variables['impactParameter'][:] - radius_of_curvature ) / 1000.0\n",
    "    d.close()\n",
    "\n",
    "    #  Plot bending angle profile. \n",
    "\n",
    "    ax.plot( bending_angle, impact_height, lw=0.2 )\n",
    "\n",
    "#  Save Figures. \n",
    "\n",
    "fname = \"all_profiles.eps\"\n",
    "print( f\"Saving figure to {fname}\" )\n",
    "fig.savefig( fname, format='eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you wish not to save the file and just display it, comment out the method \n",
    "*.savefig(...)* above. "
   ]
  },
  {
   "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}