Generic NMSIM grid results (.tig) to raster (.TIF).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convert NMSIM grid results *(.tig)* to raster *(.TIF)*\n",
    "\n",
    "We will use the library `xarray` as a means to store 4D acoustic information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ==============  import libraries  =================\n",
    "\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import interpolate\n",
    "import itertools\n",
    "import matplotlib.pyplot as plt\n",
    "import re\n",
    "import os\n",
    "import gdal\n",
    "import osr\n",
    "from tqdm.notebook import tqdm # a very helpful progress bar\n",
    "\n",
    "# ==============  define functions  =================\n",
    "\n",
    "def formatThirdOctave(raw_df):\n",
    "    \n",
    "    '''\n",
    "    Format raw chunks of SPL data from a (.tig) file into a `numpy` array; \n",
    "    we use `pandas` as a \"stepping stone\" to make this much easier.\n",
    "    '''\n",
    "    \n",
    "    # split each row of the raw data on whitespace\n",
    "    arr = np.array([t.str.split(r\"\\s+\")[0] for _, t in raw_df.iterrows()]) \n",
    "\n",
    "    # convert 2D array to dataframe, and change dtype from object to numeric\n",
    "    third = pd.DataFrame(arr).apply(pd.to_numeric, errors='coerce')  \n",
    "\n",
    "    # use the sequence column as an index instead\n",
    "    third = third.set_index(third.columns[0], drop=True) \n",
    "\n",
    "    # remove the last column (which comes from the newline character at the end of each row)\n",
    "    third = third.drop(third.columns[-1], axis=1) \n",
    "\n",
    "    # where applicable, convert from centibels to decibels (and d-prime to decimal)\n",
    "    third.loc[:, 2:] = third.loc[:, 2:]/10 \n",
    "\n",
    "    third[third == -99.9] = np.nan # use NaN for no data cells\n",
    "\n",
    "    return third.values\n",
    "\n",
    "\n",
    "def LAx(results, x=10):\n",
    "    \n",
    "    '''\n",
    "    Return exceedance level LAx,* i.e., \"the level exceeded x% of the time\",\n",
    "    and the \"*\"\" indicates an undefined integration time (in our case, it should\n",
    "    usually be about 1-2 seconds.)\n",
    "    '''\n",
    "    \n",
    "    # convert exceedance to quantile\n",
    "    q = (100 - x)/100\n",
    "    \n",
    "    values = results.loc[:,:,'LAeq'].astype('float').quantile(q, dim=\"timeStep\")\n",
    "    \n",
    "    # it's necessary to replace np.nan for 'cubic' interpolation\n",
    "    # unfortunately this affects the interpolation results (minorly)\n",
    "    # so we'll use the minimum value already existing in the array \n",
    "    min_existing = values[~np.isnan(values)].min().values\n",
    "    values = np.nan_to_num(values, nan=min_existing) \n",
    "    \n",
    "    percentile = \"{0:.0f}th percentile\".format(100 - x)\n",
    "    \n",
    "    # also include a few metadata on this metric\n",
    "    meta = {'full': percentile+\" $LA_{eq,*}$\",\n",
    "            'alias': \"LA\"+\"{0:.0f}\".format(x),\n",
    "            'unit': \"dB\"}\n",
    "    \n",
    "    return values, meta\n",
    "\n",
    "\n",
    "def LTx(results, x=10):\n",
    "    \n",
    "    '''\n",
    "    Return exceedance level LTx,* i.e., \"the level exceeded x% of the time\",\n",
    "    and the \"*\" indicates an undefined integration time (in our case, it should\n",
    "    be close to one second. Or... then again, maybe two.)\n",
    "    \n",
    "    The 'T' stands for 'truncated' or 'traffic'; it is the band from 12.5 - 1250 Hz.\n",
    "    \n",
    "    '''\n",
    "\n",
    "    # select levels in the truncated (or 'traffic') bands; convert to [relative] pressures\n",
    "    pressures = np.power(10, results.loc[:,:,\"12.5\":\"1250\"].astype('float')/10)\n",
    "     \n",
    "    # sum across bands; convert back to levels\n",
    "    levels = 10*np.log10(pressures.sum(dim=\"thirdOct\"))\n",
    "    \n",
    "    # lots of -inf; set to a very small pressure, instead\n",
    "    levels = np.nan_to_num(levels, neginf=-99.9) \n",
    "    \n",
    "    # convert exceedance to quantile\n",
    "    q = 100 - x\n",
    "    \n",
    "    # we already have the truncated sum; now take the appropriate quantile in the time dimension\n",
    "    values = np.nanpercentile(levels.astype('float'), q, axis=1)\n",
    "    \n",
    "    # it's necessary to replace np.nan for 'cubic' interpolation\n",
    "    # unfortunately this affects the interpolation results (minorly)\n",
    "    # so we'll use the minimum value already existing in the array \n",
    "    min_existing = values[~np.isnan(values)].min()\n",
    "    values = np.nan_to_num(values, nan=min_existing, neginf=min_existing)\n",
    "    \n",
    "    percentile = \"{0:.0f}th percentile\".format(100 - x)\n",
    "    \n",
    "    # also include a few metadata on this metric\n",
    "    meta = {'full': percentile+\" $LT_{eq,*}$\",\n",
    "            'alias': \"LT\"+\"{0:.0f}\".format(x),\n",
    "            'unit': \"dB\"}\n",
    "    \n",
    "    return values, meta\n",
    "\n",
    "\n",
    "def TimeAbove(results, threshold=35.0, return_as_time=True):\n",
    "    \n",
    "    '''\n",
    "    Tabulate the time equal-to or exceeding a sound level threshold for each grid reciever point.\n",
    "    \n",
    "    \n",
    "    Inputs\n",
    "    ------\n",
    "    results: (xarray object) 3D DataArray containing .tig results\n",
    "    \n",
    "    threshold: (float) a sound level threshold\n",
    "    \n",
    "    return_as_time: (bool) if True - return the absolute amount of time above the threshold, in seconds; \n",
    "                           if False - return the percentage of the model duration above the threshold\n",
    "                           \n",
    "    Returns\n",
    "    -------\n",
    "    [if `return_as_time` is True] TimeAbove: (numpy array) duration above threshold for each grid cell (in seconds)\n",
    "    [if `return_as_time` is False] pTimeAbove: (numpy array) percentage of total duration \n",
    "                                                             above threshold for each grid cell\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    # we'll work with LAeq, here\n",
    "    values = results.loc[:,:,'LAeq'].astype('float')\n",
    "    \n",
    "    # each time step during the model is slightly different (annoying but true)\n",
    "    timeStepDurations = np.diff(results.loc[:,:,\"time_s\"])\n",
    "    \n",
    "    # total duration for each grid site (should be nearly identical)\n",
    "    timeStepTotals = np.sum(timeStepDurations, axis=1)\n",
    "    \n",
    "    # for each cell in the matrix, return a boolean mask of threshold exceedance\n",
    "    above_bools = (values >= threshold)[:, :-1]\n",
    "    \n",
    "    # the duration (in seconds) that the threshold was exceeded during the model run\n",
    "    TimeAbove = np.array([np.sum(d[b]) for d, b in zip(timeStepDurations, above_bools)]) \n",
    "        \n",
    "    if(return_as_time):\n",
    "        \n",
    "        # also include a few metadata on this metric\n",
    "        meta = {'full': \"Time Above $LA_{eq,*} =$\" + \" {0:.1f} dB\".format(threshold),\n",
    "                'alias': \"TA\"+\"{0:.0f}\".format(threshold),\n",
    "                'unit': \"Time Above (seconds)\"}\n",
    "\n",
    "        return TimeAbove, meta\n",
    "\n",
    "    else:\n",
    "        \n",
    "        # also include a few metadata on this metric\n",
    "        meta = {'full': \"Time Above $LA_{eq,*} =$\" + \" {0:.1f} dB\".format(threshold),\n",
    "                'alias': \"TA\"+\"{0:.0f}_percent\".format(threshold),\n",
    "                'unit': \"Time Above (% of event duration)\"}\n",
    "        \n",
    "        # percentage: time above over total time\n",
    "        pTimeAbove = 100*TimeAbove/timeStepTotals\n",
    "        \n",
    "        return pTimeAbove, meta\n",
    "\n",
    "    \n",
    "def SEL(results):\n",
    "    \n",
    "    '''\n",
    "    Tabulate the Sound Exposure Level (SEL) for each grid cell. SEL is a way of measuring \n",
    "    the overall 'noise dose' at a location.\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    # we'll work with LAeq\n",
    "    # ...and, because we use time differencing, we need to drop the last value\n",
    "    values = results.loc[:,:,'LAeq'].astype('float')[:,:-1]\n",
    "\n",
    "    # each time step during the model is slightly different (annoying but true)\n",
    "    timeStepDurations = np.diff(results.loc[:,:,\"time_s\"])\n",
    "    \n",
    "    # compute the Sound Exposure Level (SEL) for each point\n",
    "    SEL = 10*np.log10(np.sum(timeStepDurations*np.power(10, values/10), axis=1).astype('float'))\n",
    "\n",
    "    # it's necessary to replace -inf for 'cubic' interpolation\n",
    "    # unfortunately this affects the interpolation results (minorly)\n",
    "    # so we'll use the minimum value already existing in the array \n",
    "    min_existing = SEL[~np.isnan(SEL) & np.isfinite(SEL)].min()\n",
    "    SEL = np.nan_to_num(SEL, nan=min_existing, neginf=min_existing)\n",
    "    \n",
    "    # also include a few metadata on this metric\n",
    "    meta = {'full': \"Sound Exposure Level\",\n",
    "            'alias': \"SEL\",\n",
    "            'unit': \"dB\"}\n",
    "\n",
    "    return SEL, meta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Step 1] where is the NMSIM grid result *(.tig)* you'd like to process?\n",
    "\n",
    "#### <font color=salmon size=4>You must edit this cell.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# tig = r\"V:\\NMSIM_Alpha_Tests_AH\\NMSIM_to_TRJ_Test\\ROMOTEST\\Output_Data\\TIG_TIS\\ROMOTEST_C207_1200ft_AGL_120X120.tig\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Step 2] Pre-digest the (.tig) file to determine the data bounds\n",
    "\n",
    "This is relies heavily on Damon Joyce's R function `ConvertTIG2RDATA`, lines 18 - 59. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(tig) as f:\n",
    "\n",
    "    # we'll iterate through each line...\n",
    "    sites = []  # ...gleaning the UTM coordinates\n",
    "    site_rows = np.array([]) # ...gleaning the location of every site header\n",
    "    \n",
    "    for i, line in enumerate(f):\n",
    "\n",
    "        if(re.findall(\"---End File Header---\", line)):\n",
    "\n",
    "            # the file header only occurs once; we want to know where this block ends\n",
    "            fileHeaderLen = i\n",
    "\n",
    "        elif(re.findall(\"Site:\", line)):\n",
    "\n",
    "            # site rows occur repeatedly; document each occurrence\n",
    "            site_rows = np.append(site_rows, i)\n",
    "\n",
    "        elif(re.findall(r\"^UTM\", line)):\n",
    "\n",
    "            # glean the UTM zone that NMSIM used\n",
    "            UTM_zone = int(line.split(r\"  \")[1])\n",
    "\n",
    "            # glean the UTM coordinates for the site \n",
    "            # and append them to the list as an (x, y) pair\n",
    "            sites.append([int(line[27:33]), int(line[48:55])])\n",
    "\n",
    "\n",
    "    # convert the list-of-lists to a numpy array\n",
    "    sites = np.array(sites)\n",
    "\n",
    "# how big is each data block?\n",
    "chunk_len = np.diff(site_rows)[0] - 10\n",
    "\n",
    "# tabulate all header rows; when loading the whole file we'll tell pandas to skip these\n",
    "fileHeaderRows = np.arange(0, fileHeaderLen + 1, dtype=\"int\")\n",
    "siteHeaderRows = np.array([np.arange(s - 1, s + 9, 1) for s in site_rows], dtype=\"int\")\n",
    "\n",
    "# these are all the non-SPL rows\n",
    "rowsToSkip = np.append(fileHeaderRows, siteHeaderRows)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### An (optional) diagnostic plot showing each 'reciver site'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(sites.T[0], sites.T[1], ls=\"\", marker=\"o\", color=\"k\", ms=1)\n",
    "plt.title(\"reciever grid site locations (UTM)\", loc=\"left\", y=1.03)\n",
    "\n",
    "# for some reason matplotlib insists on converting large values to scientific notation\n",
    "# this supresses default behavior\n",
    "ax = plt.gca()\n",
    "ax.ticklabel_format(style=\"plain\")\n",
    "\n",
    "plt.xlabel(\"longitude (m)\", labelpad=15)\n",
    "plt.ylabel(\"latitude (m)\", labelpad=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Step 3] Read the data and organize into a 3D `xarray` object 😊 \n",
    "\n",
    "Technically the original data are **4D**: an $(x, y)$ spatial matrix of $(t, f)$ pressure matrices. <br><br>\n",
    "We map two spatial dimensions to a single site index $(x, y) \\rightarrow j$. Axis label:  `site`. <br>\n",
    "Because the exact timing at each point is subtly different, the second axis is labelled: `timeStep`. <br>\n",
    "The third axis is frequency in one-third octave bands, labelled: `thirdOct`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Read SPL data from the .tig as chunks\n",
    "# (...unfortunately they'll come in as a single column. UGLY.)\n",
    "dfs = pd.read_csv(tig, \n",
    "                  header=None,\n",
    "                  skiprows=rowsToSkip,\n",
    "                  skipinitialspace=True,\n",
    "                  chunksize=chunk_len)\n",
    "\n",
    "# So we format SPL data for each site into a 3D numpy array\n",
    "# all values in the array must be of the same dtype, so we'll use \"object\"\n",
    "# (It's S L O W: a `tqdm` progress bar is implemented as a courtesy...)\n",
    "SPL_data = np.array([formatThirdOctave(chunk) for chunk in tqdm(dfs, total=len(sites))][:-1],\n",
    "                    dtype=\"object\")\n",
    "\n",
    "# dim_0 corresponds to each site number [label: \"site\"]\n",
    "items = np.arange(1, SPL_data.shape[0]+1)\n",
    "\n",
    "# dim_1 (major axis) corresponds to the time steps [label: \"timeStep\"]\n",
    "major_axis = np.arange(1, SPL_data.shape[1]+1)\n",
    "\n",
    "# dim_2 (minor_axis) corresponds to the following columns [label: \"thirdOct\"]\n",
    "minor_axis = [\"time_s\", \"Leq\", \"LAeq\", \"10\", \"12.5\", \"15.8\", \"20\", \"25\", \"31.5\",\n",
    "                   \"40\", \"50\", \"63\", \"80\", \"100\", \"125\", \"160\", \"200\", \"250\", \"315\",\n",
    "                   \"400\", \"500\", \"630\", \"800\", \"1000\", \"1250\", \"1600\", \"2000\", \"2500\",\n",
    "                   \"3150\", \"4000\", \"5000\", \"6300\", \"8000\", \"10000\", \"12500\", \"d'\"]\n",
    "\n",
    "# create the xarray object\n",
    "results = xr.DataArray(SPL_data, [items, major_axis, minor_axis], \n",
    "                       dims=[\"site\",\"timeStep\",\"thirdOct\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Step 4] Compute acoustic metrics from the `xarray` object\n",
    "<font size=3>These will result in a **2D** matrix $(x, y)$ representing the metric at each reciever point.</font>\n",
    "\n",
    "*(Again, this often relies on Damon's original work in R.)*\n",
    "\n",
    "#### <font color=salmon size=4>You must edit this cell.</font>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values, meta = LAx(results, x=50)\n",
    "# values, meta = LTx(results, x=90)\n",
    "# values, meta = TimeAbove(results, threshold=35.0, return_as_time=True)\n",
    "# values, meta = SEL(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Step 5] save results as a raster *(.TIF)* using `gdal`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_size = np.percentile(np.diff(sites.T[1]), 50) # the typical difference in latitude\n",
    "xmin, xmax = sites.T[0].min(), sites.T[0].max()\n",
    "ymin, ymax = sites.T[1].min(), sites.T[1].max()\n",
    "\n",
    "# here we compute the x and y steps necessary to make a grid of square cells\n",
    "# (interpolation finer than 'site' spacing of the .tig is ill-advised; we'll approximate)\n",
    "xsteps = np.ceil((xmax - xmin)/cell_size).astype('int')\n",
    "ysteps = np.ceil((ymax - ymin)/cell_size).astype('int')\n",
    "\n",
    "# then use the function `mgrid` to create new pairs of points \n",
    "# at each pair we will interpolate the expected value\n",
    "xy_new = np.mgrid[xmin:xmax:complex(0, xsteps), \n",
    "                  ymin:ymax:complex(0, ysteps)].reshape(2,-1).T\n",
    "\n",
    "# update the cell size for when we save out the raster (it'll be slightly different than the .tig)\n",
    "cell_size_final = np.percentile([np.linalg.norm(xy_new[i+1] - xy_new[i]) for \n",
    "                                              i in np.arange(len(xy_new)-1)], 50)\n",
    "\n",
    "# we can interpolate new values on a grid\n",
    "# NOTE: 'cubic' would probably provide finer results, but the nan values cause problems\n",
    "interpolated = interpolate.griddata(sites, values, xy_new, method='cubic')\n",
    "\n",
    "# reshape interpolated values into a raster and re-orient\n",
    "raster = np.reshape(interpolated, (xsteps, ysteps))\n",
    "raster = np.flip(np.rot90(raster, k=3), axis=1)\n",
    "\n",
    "# a preview of the results! handy\n",
    "plt.figure(figsize=(8, 7))\n",
    "plt.imshow(raster, origin=\"lower\", aspect=\"auto\", cmap=\"magma\")\n",
    "plt.suptitle(meta['full'], ha=\"left\", fontsize=16, x=0.12, y=0.96)\n",
    "plt.title(os.path.basename(tig), loc=\"left\")\n",
    "cb = plt.colorbar(shrink=0.8)\n",
    "cb.set_label(meta['unit'], labelpad=15, fontsize=13)\n",
    "# plt.savefig(os.path.expanduser(\"~/Desktop\") + os.sep + meta['alias'] + \"_.png\",\n",
    "#             dpi=150, bbox_inches=\"tight\")\n",
    "plt.show()\n",
    "\n",
    "# epsg codes for UTM zones in the United States\n",
    "epsg_lookup = {1:26901, 2:26902, 3:26903, 4:26904, 5:26905, 6:26906, \n",
    "               7:26907, 8:26908, 9:26909, 10:26910, 11:26911, 12:26912, \n",
    "               13:26913, 14:26914, 15:26915, 16:26916, 17:26917,\n",
    "               18:26918, 19:26919, 20:26920}\n",
    "\n",
    "# initialize a spatial reference and look up the appropriate EPSG code as integer\n",
    "srs = osr.SpatialReference()\n",
    "srs.ImportFromEPSG(int(epsg_lookup[UTM_zone]))\n",
    "\n",
    "# path to the output raster\n",
    "out_path = tig[:-4] + \"_\" + meta['alias'] + \".TIF\"\n",
    "\n",
    "# write the raster to an array\n",
    "driver = gdal.GetDriverByName('GTiff')\n",
    "out = driver.Create(out_path, raster.shape[1], raster.shape[0], 1, gdal.GDT_Float32)\n",
    "out.GetRasterBand(1).WriteArray(raster)\n",
    "\n",
    "out.SetGeoTransform((xmin, cell_size_final, 0, ymin, 0, cell_size_final))\n",
    "out.SetProjection(srs.ExportToWkt())\n",
    "out.FlushCache()\n",
    "out=None"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}