{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vmlzyjyZStZ3"
   },
   "source": [
    "   ![alt text](https://www.mbari.org/wp-content/uploads/2014/11/logo-mbari-3b.png \"MBARI\")\n",
    "\n",
    "  <div align=\"left\">Copyright (c) 2021, MBARI</div>\n",
    "    \n",
    "  * Distributed under the terms of the GPL License\n",
    "  * Maintainer: dcline@mbari.org\n",
    "  * Authors: Danelle Cline dcline@mbari.org, John Ryan ryjo@mbari.org"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6h0JRTviT-fW"
   },
   "source": [
    "## Basic Exploration of the 16 kHz Pacific Ocean Audio Data in the AWS Open Data Registry\n",
    "\n",
    "---\n",
    "An extensive (5+ years and growing) archive of sound recordings from a deep-sea location [along the eastern margin of the North Pacific Ocean](https://www.mbari.org/at-sea/cabled-observatory/) has been made available through AWS Open data.  Temporal coverage of the recording archive has been 95% since project inception in July 2015.  The original recordings have a sample rate of 256 kHz.  For many research applications it is convenient to work with data having a lower sample rate.  This notebook illustrates basic methods to access and process a calibrated spectrogram from the decimated 16 kHz audio archive.\n",
    "\n",
    "If you use this data set, please **[cite our project](https://ieeexplore.ieee.org/document/7761363).**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "FpDVTEOEZ-qn"
   },
   "source": [
    "## Data Overview\n",
    "The decimated audio data are in [WAV](https://en.wikipedia.org/wiki/WAV) format in an s3 bucket named <b>pacific-sound-16khz</b>.  They are further organized by year and month.  Buckets are stored as objects, so the data isn't physically stored in folders or directories as you may be famaliar with, but you can think of it conceptually as follows:\n",
    "\n",
    "```\n",
    "pacific-sound-16khz\n",
    "      |\n",
    "      ----2020\n",
    "        |\n",
    "        |----01\n",
    "        ...\n",
    "        |----12\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "H8P2AGrQUOxu"
   },
   "source": [
    "## Install dependencies\n",
    "First, let's install dependencies and include all packages used in this tutorial. This only needs to be done once for the duration of this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bu5cv0WFx3hX"
   },
   "outputs": [],
   "source": [
    "!pip install -q boto3\n",
    "!pip install -q soundfile\n",
    "!pip install -q scipy\n",
    "!pip install -q numpy\n",
    "!pip install -q matplotlib\n",
    "import boto3, botocore\n",
    "from botocore import UNSIGNED\n",
    "from botocore.client import Config\n",
    "from six.moves.urllib.request import urlopen\n",
    "from pathlib import Path\n",
    "import io\n",
    "import os\n",
    "import scipy\n",
    "from scipy import signal\n",
    "import numpy as np\n",
    "import soundfile as sf\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "87EGQ0Oc5oD5"
   },
   "source": [
    "## List the contents of a monthly directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kgsKVYlPduHI"
   },
   "outputs": [],
   "source": [
    "s3_client = boto3.client('s3',\n",
    "    aws_access_key_id='',\n",
    "    aws_secret_access_key='', \n",
    "    config=Config(signature_version=UNSIGNED))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "_t5aG6C66KZD"
   },
   "outputs": [],
   "source": [
    "year = \"2018\"\n",
    "month = \"01\"\n",
    "bucket = 'pacific-sound-16khz'\n",
    "\n",
    "for obj in s3_client.list_objects_v2(Bucket=bucket, Prefix=f'{year}/{month}')['Contents']:\n",
    "    print(obj['Key'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "636LHKTbAgNf"
   },
   "source": [
    "## Retrieve metadata for a file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "X2y_ADUFWOKP"
   },
   "outputs": [],
   "source": [
    "year = \"2018\"\n",
    "month = \"01\"\n",
    "filename = 'MARS-20180101T000000Z-16kHz.wav'\n",
    "bucket = 'pacific-sound-16khz'\n",
    "key = f'{year}/{month}/{filename}'\n",
    "\n",
    "url = f'https://{bucket}.s3.amazonaws.com/{key}'\n",
    "\n",
    "sf.info(io.BytesIO(urlopen(url).read(20_000)), verbose=True) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dHtzJHKtuCti"
   },
   "source": [
    "## Calibrated Spectrum Levels\n",
    "The file metadata retrieved above includes a link to a json file containing complete deployment metadata.  For convenience, those links (one for each hydrophone deployment) are: \n",
    "*   https://bitbucket.org/mbari/pacific-sound/src/master/MBARI_MARS_Hydrophone_Deployment01.json\n",
    "*   https://bitbucket.org/mbari/pacific-sound/src/master/MBARI_MARS_Hydrophone_Deployment02.json\n",
    "\n",
    "Frequency-dependent hydrophone sensitivity data are included in these json files.  Considering the 16 kHz data, hydrophone sensitivity at only the two lowest frequencies are relevant.  For the example file identified above, which is from the second hydrophone deployment, these sensitivities at 250 Hz and 10 kHz are -177.90 and -176.80 dB V re $\\mu$Pa, respectively.  For simplicity in this example, we can apply the average sensitivity, -177.35 dB V re $\\mu$Pa.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WGHEshKmDGI9"
   },
   "source": [
    "### Save the daily file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Xox8b2IDh5PZ"
   },
   "outputs": [],
   "source": [
    "# The 16 kHz daily file may be too large to fit into memory for some instances, so let's download it first\n",
    "# then seek into the file to read a chunk later\n",
    "print(f'Copying s3://{bucket}/{key}')\n",
    "  \n",
    "# only download if needed\n",
    "if not Path(filename).exists():\n",
    "    s3 = boto3.resource('s3',\n",
    "        aws_access_key_id='',\n",
    "        aws_secret_access_key='',\n",
    "        config=Config(signature_version=UNSIGNED))\n",
    "\n",
    "    # Alternatively, it can be downloaded directly in SageMaker with\n",
    "    # !aws s3 cp s3://{bucket}/{key} . \n",
    "\n",
    "    print('Downloading') \n",
    "    s3.Bucket(bucket).download_file(key, filename)\n",
    "    print('Done')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read hours 9-12 of the day"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_rate = int(16e3)\n",
    "start_frame = int(sample_rate*9*3600)\n",
    "duration_frames =  int(sample_rate*3*3600)\n",
    "\n",
    "pacsound_file = sf.SoundFile(filename)\n",
    "pacsound_file.seek(start_frame)\n",
    "x = pacsound_file.read(duration_frames, dtype='float32')    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "LKivxTATjYWL"
   },
   "outputs": [],
   "source": [
    "v = x*3   # convert scaled voltage to volts\n",
    "v.shape, v.size, sample_rate\n",
    "a = np.arange(v.size)+1\n",
    "# define segment processing\n",
    "nsec = (v.size)/sample_rate # number of seconds in vector\n",
    "spa = 60  # seconds per average\n",
    "nseg = int(nsec/spa)\n",
    "print(nseg,'segments of length',spa,'seconds','in',nsec,'seconds of audio')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Sn0ez0DzkPP-"
   },
   "outputs": [],
   "source": [
    "# initialize empty spectrogram matrix\n",
    "nfreq = int(sample_rate/2+1)\n",
    "nfreq,nseg\n",
    "sg = np.empty((nfreq, nseg), float)\n",
    "sg.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "eu2i-CF4kWTV"
   },
   "outputs": [],
   "source": [
    "# get window for welch\n",
    "w = scipy.signal.get_window('hann',sample_rate)\n",
    "\n",
    "# process spectrogram\n",
    "for x in range(0,nseg):\n",
    "  cstart = x*spa*sample_rate\n",
    "  cend = (x+1)*spa*sample_rate\n",
    "  f,psd = scipy.signal.welch(v[cstart:cend],fs=sample_rate,window=w,nfft=sample_rate)\n",
    "  psd = 10*np.log10(psd) + 177.35\n",
    "  sg[:,x] = psd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JhXzdmYONg9V"
   },
   "source": [
    "### Plot the spectrogram of the 6-hour window\n",
    "Note: The sharp drop in signal approaching 8 kHz reflects the attributes of the decimation filter applied to produce the 16 kHz data from the original 256 kHz data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NFRo06jRkiLB"
   },
   "outputs": [],
   "source": [
    "plt.figure(dpi=300)\n",
    "im = plt.imshow(sg,aspect='auto',origin='lower',vmin=30,vmax=100)\n",
    "plt.yscale('log')\n",
    "plt.ylim(10,8000)\n",
    "plt.colorbar(im)\n",
    "plt.annotate(\"humpback whale song (above ~70 Hz)\",(0,180))\n",
    "plt.annotate(\"blue whale song (B-call 3rd harmonic)\",(100,40))\n",
    "plt.annotate(\"fin whale song (20 Hz pulses)\",(100,18))\n",
    "plt.xlabel('Minute of day')\n",
    "plt.ylabel('Frequency (Hz)')\n",
    "plt.title('Calibrated spectrum levels')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5bRrM9ZcDulz"
   },
   "source": [
    "### Zoom-in using a 1-second resolution to see a more detailed of whale song\n",
    "The temporal resolution of the full-day spectrogram above is too coarse to see the structure of individual whale songs.  So, let's produce and view a calibrated spectrogram with more detail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NJlNq3yNEt6C"
   },
   "outputs": [],
   "source": [
    "# examine the 9th hour of the day at 1 second resolution\n",
    "spa = 1  # seconds per average\n",
    "nseg = 3600\n",
    "sg2 = np.empty((nfreq, nseg), float)\n",
    "start_sample = int(3600*sample_rate+1)\n",
    "sg2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "GO7iUR2FmC04"
   },
   "outputs": [],
   "source": [
    "# process spectrogram\n",
    "for x in range(0,nseg):\n",
    "  cstart = start_sample + x*spa*sample_rate\n",
    "  cend = cstart+spa*sample_rate\n",
    "  f,psd = scipy.signal.welch(v[cstart:cend],fs=sample_rate,window=w,nfft=sample_rate)\n",
    "  psd = 10*np.log10(psd) + 177.35\n",
    "  sg2[:,x] = psd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "TgmlT_-nGyR2"
   },
   "outputs": [],
   "source": [
    "# plot it\n",
    "plt.figure(dpi=300)\n",
    "im = plt.imshow(sg2,aspect='auto',origin='lower',vmin=30,vmax=100)\n",
    "plt.yscale('log')\n",
    "plt.ylim(10,8000)\n",
    "plt.colorbar(im)\n",
    "plt.annotate(\"humpback whale song\",(0,180))\n",
    "plt.annotate(\"blue whale song\",(1200,40))\n",
    "plt.annotate(\"fin whale song\",(1400,18))\n",
    "plt.annotate(\"earthquake\",(150,11))\n",
    "plt.xlabel('Second of hour')\n",
    "plt.ylabel('Frequency (Hz)')\n",
    "plt.title('Calibrated spectrum levels')"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "machine_shape": "hm",
   "name": "PacificSound_16kHz.ipynb",
   "private_outputs": true,
   "provenance": []
  },
  "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.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}