{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt, ticker\n",
    "from scipy import stats\n",
    "from mintpy.objects import timeseries, sensor\n",
    "from mintpy.utils import readfile, utils as ut, plot as pp\n",
    "\n",
    "\n",
    "##---------- plot range delay time series\n",
    "def adjust_ts_axis_format(ax, x, vstep=0.3, xlim=None, ylim=None, ylabel='cumulative\\nslant range offset [m]', xticklabels=True):\n",
    "    ax.tick_params(which='both', direction='in', bottom=True, top=True, left=True, right=True)\n",
    "    ax.grid('on', linestyle='--')\n",
    "    ax.yaxis.set_major_locator(ticker.MultipleLocator(vstep))\n",
    "    pp.auto_adjust_xaxis_date(ax, x, every_year=1, buffer_year=None)\n",
    "    if xlim is not None:\n",
    "        ax.set_xlim(xlim)\n",
    "    if ylim is not None:\n",
    "        ax.set_ylim(ylim)\n",
    "    if ylabel:\n",
    "        ax.set_ylabel(ylabel)\n",
    "    if not xticklabels:\n",
    "        ax.set_xticklabels([])\n",
    "    return\n",
    "\n",
    "def plot_dot_figure(ax, x, y, vlim, vstep=0.2, xname=None, yname=None, fbase='offset'):\n",
    "    # omit nan value\n",
    "    flag = (~np.isnan(x)) * (~np.isnan(y))\n",
    "    x = x[flag]\n",
    "    y = y[flag]\n",
    "\n",
    "    ax.plot(x, y, 'k.')\n",
    "    ax.plot(vlim, vlim, 'k--', lw=1)\n",
    "\n",
    "    # axis format\n",
    "    ax.tick_params(which='both', direction='in', bottom=True, top=True, left=True, right=True)\n",
    "    ax.xaxis.set_major_locator(ticker.MultipleLocator(vstep))\n",
    "    ax.yaxis.set_major_locator(ticker.MultipleLocator(vstep))\n",
    "    #ax.set_yticklabels([])\n",
    "    ax.yaxis.tick_right()\n",
    "    ax.yaxis.set_label_position(\"right\")\n",
    "    ax.set_xlim(vlim)\n",
    "    ax.set_ylim(vlim)\n",
    "    #ax.set_xlabel('${{{b}}}_{{{n}}}$ [m]'.format(b=fbase, n=xname), color='C0')\n",
    "    #ax.set_ylabel('${{{b}}}_{{{n}}}$ [m]'.format(b=fbase, n=yname), color='C1', labelpad=ylabelpad)\n",
    "    ax.set_aspect('equal', 'box')\n",
    "\n",
    "    # stats\n",
    "    rmse = np.sqrt(np.sum((y - x)**2) / (x.size - 1))\n",
    "    r2 = stats.linregress(y, x)[2]\n",
    "    msg  = r'$R^2$ = {:.2f}'.format(r2) + '\\n'\n",
    "    msg += r'$RMSE$ = {:.1f} cm'.format(rmse*100) + '\\n'\n",
    "    #ax.annotate(msg, xy=(0.05, 0.65), xycoords='axes fraction', color='k', ha='left')\n",
    "    #ax.annotate(msg, xy=(0.95, 0.07), xycoords='axes fraction', color='k', ha='right')\n",
    "    print(msg)\n",
    "\n",
    "    return\n",
    "\n",
    "##---------- read list of time series files for RMSE analysis\n",
    "def read_ts_files(fnames, print_msg=True, print_max=True, rm_SenD_outlier=True):\n",
    "    fnames = [x for x in fnames if os.path.isfile(x)]\n",
    "    proj_dir = os.path.dirname(fnames[0])\n",
    "    proj_name = os.path.basename(os.path.dirname(proj_dir))\n",
    "\n",
    "    ### geometry files\n",
    "    geom_file = os.path.join(proj_dir, 'inputs/geometryRadar.h5')\n",
    "    mask_file = os.path.join(proj_dir, 'maskResInv.h5')\n",
    "\n",
    "    # date info\n",
    "    date_list = timeseries(fnames[0]).get_date_list()\n",
    "    num_date = len(date_list)\n",
    "\n",
    "    # for ChileSenDT156, remove the outlier at 2015-04-02\n",
    "    flag = np.ones(num_date, dtype=np.bool_)\n",
    "    if rm_SenD_outlier and proj_name == 'ChileSenDT156':\n",
    "        flag[date_list.index('20150402')] = 0\n",
    "\n",
    "    # lalo --> box\n",
    "    lalo = [-21.30, -67.39] if proj_name.startswith('ChileSen') else None\n",
    "    if lalo is None:\n",
    "        box = None\n",
    "    else:\n",
    "        atr = readfile.read_attribute(fnames[0])\n",
    "        coord = ut.coordinate(atr, lookup_file=geom_file)\n",
    "        y, x = coord.geo2radar(lalo[0], lalo[1])[:2]\n",
    "        win = 10\n",
    "        box = (x-win, y-win, x+win+1, y+win+1)\n",
    "\n",
    "    # read data into lists\n",
    "    tsDict, rDict = dict(), dict()\n",
    "    mask = readfile.read(mask_file, box=box)[0].flatten()\n",
    "    for i, fname in enumerate(fnames):\n",
    "        # label\n",
    "        fbase = os.path.basename(fname).replace('timeseriesRg', 'SAR').replace('.h5', '')\n",
    "        label = ' - '.join(fbase.split('_'))\n",
    "\n",
    "        # read data\n",
    "        ts_data = readfile.read(fname, box=box)[0].reshape(num_date, -1)[flag][:, mask]\n",
    "        ts_med = np.nanmedian(ts_data, axis=-1)\n",
    "        ts_med -= np.nanmedian(ts_med)\n",
    "\n",
    "        # print out stats\n",
    "        rmse = ut.root_mean_sq_error(ts_med) * 100\n",
    "        dmax = np.nanmax(np.abs(ts_med)) * 100\n",
    "        if print_msg:\n",
    "            msg = ''\n",
    "            if i == 0:\n",
    "                msg += f'{proj_name}: RMSE / MAX\\n' if print_max else f'{proj_name}: RMSE\\n'\n",
    "            #fdigit = 0 if 'Alos2' in proj_name else 1\n",
    "            fdigit = 1\n",
    "            msg += '    {l:<40}: {r:6.{f}f}'.format(l=label, r=rmse, f=fdigit)\n",
    "            msg += ' / {:4.0f}'.format(dmax) if print_max else ''\n",
    "            msg += ' cm'\n",
    "            print(msg)\n",
    "\n",
    "        # save for outputs\n",
    "        tsDict[label] = ts_med\n",
    "        rDict[label] = rmse\n",
    "\n",
    "    return proj_name, tsDict, rDict"
   ]
  },
  {
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}