SeismoSocialDistancing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ground Motion Displacement RMS vs Time\n",
    "\n",
    "*an example simple tutorial for getting seismic data, computing the power spectral densities, extracting the RMS and plotting*\n",
    "\n",
    "Required:\n",
    "\n",
    "- python\n",
    "- obspy (and its dependencies)\n",
    "- pandas\n",
    "- jupyter\n",
    "- notebook\n",
    "- tqdm\n",
    "\n",
    "this should be easy to set up in a conda env: ``conda create -c conda-forge -n covid python=3.7 obspy pandas jupyter notebook tqdm``\n",
    "\n",
    "Author: Thomas Lecocq @seismotom, Fred Massin @fmassin, Claudio Satriano @claudiodsf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:46:13.572490Z",
     "start_time": "2020-05-08T06:45:59.783138Z"
    }
   },
   "outputs": [],
   "source": [
    "import datetime\n",
    "import os\n",
    "from glob import glob\n",
    "\n",
    "import matplotlib\n",
    "matplotlib.rcParams['pdf.fonttype'] = 42  # to edit text in Illustrator\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "import matplotlib.patheffects as pe\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import tqdm\n",
    "import warnings\n",
    "\n",
    "from obspy import UTCDateTime, read\n",
    "from obspy.clients.fdsn import Client\n",
    "from obspy.clients.fdsn.client import FDSNNoDataException\n",
    "from obspy.signal import PPSD\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import seismosocialdistancing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Define Start/End dates and Seismic Channel\n",
    "\n",
    "You'll have to make sure the seed_id you request is indeed available from the ``data_provider``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:46:14.089189Z",
     "start_time": "2020-05-08T06:46:13.580338Z"
    }
   },
   "outputs": [],
   "source": [
    "# Make sure you take at least a full week (>=7 days) before the first \"ban\"\n",
    "start = UTCDateTime(\"2020-03-01\")\n",
    "# Leaving UTCDateTime() empty means \"now\":\n",
    "end = UTCDateTime(\"2020-05-01\")\n",
    "\n",
    "network = \"BE\"\n",
    "station = \"UCC\"\n",
    "location = \"\"\n",
    "channel = \"HHZ\"\n",
    "dataset = \"example\"\n",
    "time_zone = \"Europe/Brussels\"\n",
    "sitedesc = \"in Uccle (Brussels, BE)\"\n",
    "\n",
    "data_provider = \"ODC\"\n",
    "logo = None # 'https://upload.wikimedia.org/wikipedia/commons/thumb/4/44/Logo_SED_2014.png/220px-Logo_SED_2014.png'\n",
    "bans = {\"2020-03-15 00:00\":'Restaurants/Bars/Schools closed', \n",
    "        \"2020-03-18 12:00\":'Non-essential shops closed'}\n",
    "\n",
    "datelist = pd.date_range(start.datetime, min(end, UTCDateTime()).datetime, freq=\"D\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Download the seismic waveform data\n",
    "\n",
    "This step is coded so that only the last day is redownloaded if the daily files are present on the disk.\n",
    "\n",
    "The request gets the target day +- 30 minutes to avoid having gaps at the end of each day (need 1 window covering midnight)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:46:17.701310Z",
     "start_time": "2020-05-08T06:46:14.120096Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "c = Client(data_provider)\n",
    "\n",
    "nslc = \"{}.{}.{}.{}\".format(network, station, location, channel)\n",
    "# make sure that wildcard characters are not in nslc\n",
    "nslc = nslc.replace(\"*\", \"\").replace(\"?\", \"\")\n",
    "pbar = tqdm.tqdm(datelist)\n",
    "for day in pbar:\n",
    "    datestr = day.strftime(\"%Y-%m-%d\")\n",
    "    fn = \"{}_{}_{}.mseed\".format(dataset, datestr, nslc)\n",
    "    if day != UTCDateTime().datetime and os.path.isfile(fn):\n",
    "        continue\n",
    "    else:\n",
    "        pbar.set_description(\"Fetching %s\" % fn)\n",
    "        try: \n",
    "            st = c.get_waveforms(network, station, location, channel,\n",
    "                                  UTCDateTime(day)-1801, UTCDateTime(day)+86400+1801,\n",
    "                                  attach_response=True)\n",
    "            # the following two lines here prevent a crash when there is no data on a day\n",
    "            if len(st) > 0:\n",
    "                st.write(fn)\n",
    "        except FDSNNoDataException:\n",
    "            pbar.set_description(\"No data on FDSN server for %s\" % fn)\n",
    "            continue\n",
    "        \n",
    "resp = c.get_stations(UTCDateTime(day), network=network, station=station, location=location,\n",
    "                      channel=channel, level=\"response\")\n",
    "print(resp)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Compute PPSDs using custom parameters\n",
    "\n",
    "These parameters are set to allow the PSDs to be \"nervous\", not as smooth as the default PQLX ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:47:10.106003Z",
     "start_time": "2020-05-08T06:46:17.724247Z"
    }
   },
   "outputs": [],
   "source": [
    "force_reprocess = False\n",
    "pbar = tqdm.tqdm(datelist)\n",
    "for day in pbar:\n",
    "    datestr = day.strftime(\"%Y-%m-%d\")\n",
    "    fn_in = \"{}_{}_{}.mseed\".format(dataset, datestr, nslc)\n",
    "    pbar.set_description(\"Processing %s\" % fn_in)\n",
    "    if not os.path.isfile(fn_in):\n",
    "        continue\n",
    "    stall = read(fn_in, headonly=True)\n",
    "    for mseedid in list(set([tr.id for tr in stall])):\n",
    "        fn_out = \"{}_{}_{}.npz\".format(dataset, datestr, mseedid)\n",
    "        if os.path.isfile(fn_out) and not force_reprocess:\n",
    "            continue\n",
    "        st = read(fn_in, sourcename=mseedid)\n",
    "        st.attach_response(resp)\n",
    "        ppsd = PPSD(st[0].stats, metadata=resp,\n",
    "                    ppsd_length=1800, overlap=0.5,\n",
    "                    period_smoothing_width_octaves=0.025,\n",
    "                    period_step_octaves=0.0125,\n",
    "                    period_limits=(0.008, 50),\n",
    "                    db_bins=(-200, 20, 0.25))\n",
    "        with warnings.catch_warnings():\n",
    "            warnings.simplefilter(\"ignore\")\n",
    "            ppsd.add(st)\n",
    "        ppsd.save_npz(fn_out[:-4])\n",
    "        del st, ppsd\n",
    "    del stall\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: Reload daily PSDs from the disk and create a single PPSD object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:47:16.462083Z",
     "start_time": "2020-05-08T06:47:10.109992Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ppsds = {}\n",
    "pbar = tqdm.tqdm(datelist)\n",
    "for day in pbar:\n",
    "    datestr = day.strftime(\"%Y-%m-%d\")\n",
    "    fn_pattern = \"{}_{}_*.npz\".format(dataset, datestr)\n",
    "    pbar.set_description(\"Reading %s\" % fn_pattern)\n",
    "    for fn in glob(fn_pattern):\n",
    "        mseedid = fn.replace(\".npz\", \"\").split(\"_\")[-1]\n",
    "        if mseedid not in ppsds:\n",
    "            ppsds[mseedid] = PPSD.load_npz(fn)#, allow_pickle=True)\n",
    "        else:\n",
    "            with warnings.catch_warnings():\n",
    "                warnings.simplefilter(\"ignore\")\n",
    "                ppsds[mseedid].add_npz(fn)#, allow_pickle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 6: Standard plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:49:46.636756Z",
     "start_time": "2020-05-08T06:47:16.467090Z"
    }
   },
   "outputs": [],
   "source": [
    "[ppsd.plot(max_percentage=10) for mseedid, ppsd in ppsds.items()]\n",
    "[ppsd.plot_temporal(0.10) for mseedid, ppsd in ppsds.items()]\n",
    "[ppsd.plot_spectrogram(clim=(-160,-100)) for mseedid, ppsd in ppsds.items()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 7: Process PSDs to extract the RMS(displacement)\n",
    "\n",
    "This can be done for multiple filters at once (``freqs`` below):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:50:06.581418Z",
     "start_time": "2020-05-08T06:49:49.777763Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Define frequency bands of interest:\n",
    "freqs = [(0.1,1.0),(1.0,20.0),(4.0,14.0),(4.0,20.0)]\n",
    "\n",
    "displacement_RMS = {}\n",
    "for mseedid, ppsd in tqdm.tqdm(ppsds.items()):\n",
    "    ind_times = pd.DatetimeIndex([d.datetime for d in ppsd.current_times_used])\n",
    "    data = pd.DataFrame(ppsd.psd_values, index=ind_times, columns=1./ppsd.period_bin_centers)\n",
    "    data = data.sort_index(axis=1)\n",
    "    displacement_RMS[mseedid] = seismosocialdistancing.df_rms(data, freqs, output=\"DISP\")\n",
    "    displacement_RMS[mseedid].to_csv(\"%s.csv\" % mseedid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Weekday / Time of day Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:50:10.366068Z",
     "start_time": "2020-05-08T06:50:06.585408Z"
    }
   },
   "outputs": [],
   "source": [
    "args = {'band':\"4.0-14.0\",       # might be None or commented (\"4.0-14.0\" per default) or any of the tupples in freqs\n",
    "        'time_zone':time_zone,   # required for clockplots\n",
    "        'sitedesc':sitedesc,     # might be None or commented\n",
    "        'logo':logo,             # might be None or commented\n",
    "        'bans':bans,             # might be None or commented\n",
    "        'save':'./',              # might be None or commented or a path \n",
    "        'unit':'nm'\n",
    "       }\n",
    "seismosocialdistancing.plot(displacement_RMS,\n",
    "                            type='timeseries',\n",
    "                            **args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:50:12.651399Z",
     "start_time": "2020-05-08T06:50:10.370011Z"
    }
   },
   "outputs": [],
   "source": [
    "seismosocialdistancing.plot(displacement_RMS,\n",
    "                            type='dailyplots',\n",
    "                            **args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:50:18.101276Z",
     "start_time": "2020-05-08T06:50:12.654353Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "seismosocialdistancing.plot(displacement_RMS,\n",
    "                            type='clockplots',\n",
    "                            **args)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Noise distribution over time of the day  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:50:22.177597Z",
     "start_time": "2020-05-08T06:50:18.104268Z"
    }
   },
   "outputs": [],
   "source": [
    "seismosocialdistancing.plot(displacement_RMS,\n",
    "                            type='clockmaps',\n",
    "                            **args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:50:26.135857Z",
     "start_time": "2020-05-08T06:50:22.180560Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "seismosocialdistancing.plot(displacement_RMS,\n",
    "                            type='gridmaps',\n",
    "                            **args)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:49:46.658012Z",
     "start_time": "2020-05-08T06:45:59.879Z"
    }
   },
   "source": [
    "## Temporary code: All your stations' colormapped plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-08T06:51:18.064939Z",
     "start_time": "2020-05-08T06:51:15.816607Z"
    }
   },
   "outputs": [],
   "source": [
    "normalize = True\n",
    "resample_freq = \"1H\"\n",
    "clip = (0.05,0.95)\n",
    "\n",
    "#concatenating in a single dataframe and converting to nm\n",
    "dRMS2D = pd.concat(displacement_RMS, axis=0) * 1e9\n",
    "for band in dRMS2D.columns:\n",
    "    g = dRMS2D.loc[:,band].unstack().T\n",
    "    g = g.clip(g.quantile(clip[0]),g.quantile(clip[1]),axis=1)\n",
    "    g = g.resample(resample_freq).median()\n",
    "    if normalize:\n",
    "        g -= g.quantile(0.01)\n",
    "        g /= g.quantile(0.99)\n",
    "        vmin, vmax = 0, 1\n",
    "    else:\n",
    "        vmin = vmax = None\n",
    "   \n",
    "    fig, ax = plt.subplots(1,1, figsize=(12,2*len(g.columns)))\n",
    "    label = '%sDisplacement (µm)' % [\"\",\"Normalized \"][normalize]\n",
    "\n",
    "    x = np.append(g.index.values, g.index.shift(1).values[-1])\n",
    "    y = np.append(g.columns, \"\")\n",
    "    plt.pcolormesh(x, y, g.T,  cmap=\"inferno\", vmin=vmin, vmax=vmax)\n",
    "    plt.colorbar(orientation='horizontal', shrink=0.3).set_label(label)\n",
    "    pos, l = plt.yticks()\n",
    "    plt.yticks(np.asarray(pos)+0.5, y)\n",
    "    plt.margins(0)\n",
    "    plt.title(\"%s Hz\" % band)\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.show()\n"
   ]
  },
  {
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}