{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Jupyter Notebook example for preparing subcatalogues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import h5py\n",
    "import pickle\n",
    "\n",
    "# Load OQ tools\n",
    "from openquake.hmtk.parsers.catalogue import CsvCatalogueParser\n",
    "from openquake.hmtk.seismicity.selector import CatalogueSelector\n",
    "from openquake.hmtk.parsers.catalogue.csv_catalogue_parser import CsvCatalogueWriter  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configuration files\n",
    "cat_pickle_filename = '~/model/catalogue/csv/catalogue.pkl'\n",
    "treg = '~/model/catalogue/classification/classified.hdf5'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "crustal\n",
      "crustal_deep\n",
      "int_prt\n",
      "slab_nht\n",
      "slab_prt\n"
     ]
    }
   ],
   "source": [
    "# Reading TR hdf5 file and creating the list of tectonic regions\n",
    "aaa = []\n",
    "f = h5py.File(treg, \"r\")\n",
    "for key in f.keys():\n",
    "    aaa.append(key)\n",
    "    alen = len(f[key])\n",
    "    print(key)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# earthquakes in the catalogue: 16553\n",
      "# earthquakes in this TR      : 10999\n",
      "Catalogue successfully written to cat_TR_crustal.csv\n",
      "# earthquakes in the catalogue: 16553\n",
      "# earthquakes in this TR      : 1212\n",
      "Catalogue successfully written to cat_TR_crustal_deep.csv\n",
      "# earthquakes in the catalogue: 16553\n",
      "# earthquakes in this TR      : 1933\n",
      "Catalogue successfully written to cat_TR_int_prt.csv\n",
      "# earthquakes in the catalogue: 16553\n",
      "# earthquakes in this TR      : 626\n",
      "Catalogue successfully written to cat_TR_slab_nht.csv\n",
      "# earthquakes in the catalogue: 16553\n",
      "# earthquakes in this TR      : 296\n",
      "Catalogue successfully written to cat_TR_slab_prt.csv\n"
     ]
    }
   ],
   "source": [
    "# for each label, create the subcatalogue\n",
    "tot_lab = np.zeros(alen)\n",
    "for label in (aaa):\n",
    "    csv_filename = \"cat_TR_%s.csv\"%(label)\n",
    "    f = h5py.File(treg,'r')\n",
    "    tr = f[label][:]\n",
    "    f.close()\n",
    "    if sum(tr) > 0:\n",
    "        tmp_lab = tr*1\n",
    "        tot_lab = tot_lab+tmp_lab\n",
    "        catalogue = pickle.load(open(cat_pickle_filename, 'rb'))\n",
    "        for lab in ['month', 'day', 'hour', 'minute', 'second']:\n",
    "            idx = np.isnan(catalogue.data[lab])\n",
    "            if lab == 'day' or lab == 'month':\n",
    "                catalogue.data[lab][idx] = 1\n",
    "            elif lab == 'second':\n",
    "                catalogue.data[lab][idx] = 0.0\n",
    "            else:\n",
    "                catalogue.data[lab][idx] = 0\n",
    "        selector = CatalogueSelector(catalogue, create_copy=False)\n",
    "        print('# earthquakes in the catalogue: {:d}'.format(len(catalogue.data['longitude'])))\n",
    "        catalogue = selector.select_catalogue(tr)\n",
    "        \n",
    "        print('# earthquakes in this TR      : {:d}'.format(len(catalogue.data['longitude'])))\n",
    "        # Sub-catalogue\n",
    "        csvcat = CsvCatalogueWriter(csv_filename) \n",
    "        # Write the purged catalogue\n",
    "        csvcat.write_file(catalogue)\n",
    "        print(\"Catalogue successfully written to %s\" % csv_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# earthquakes: 16553\n",
      "# earthquakes: 1487\n",
      "Catalogue successfully written to cat_TR_unclassified.csv\n"
     ]
    }
   ],
   "source": [
    "# also make a catalogue of unclassified earthquakes\n",
    "tr_undef = abs(tot_lab-1)\n",
    "catalogue = pickle.load(open(cat_pickle_filename, 'rb'))\n",
    "selector = CatalogueSelector(catalogue, create_copy=False)\n",
    "print('# earthquakes: {:d}'.format(len(catalogue.data['longitude'])))\n",
    "catalogue = selector.select_catalogue(tr_undef)\n",
    "print('# earthquakes: {:d}'.format(len(catalogue.data['longitude'])))\n",
    "# Sub-catalogue\n",
    "csv_filename = \"cat_TR_unclassified.csv\"\n",
    "csvcat = CsvCatalogueWriter(csv_filename) \n",
    "# Write the purged catalogue\n",
    "csvcat.write_file(catalogue)\n",
    "print(\"Catalogue successfully written to %s\" % csv_filename)"
   ]
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}