hyp3_mintpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "aa7e2621-be81-4b3d-89d7-a58ad655411a",
   "metadata": {},
   "source": [
    "# MintPy Processing\n",
    "\n",
    "This notebook performs the MintPy timeseries processing. If you you're processing an ascending and descending track you'll need to run this workflow twice. Once for each track. I've had to add some strange working directory management to get this workflow to run on my local machine (run the load_data step in the directory the hyp3 data is located in, then copy it into a processing directory) but feel free to change this workflow if your machine is different."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17c93a60-b988-450b-a627-9c61df4305e3",
   "metadata": {
    "tags": []
   },
   "source": [
    "## setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61a13a80-084f-45be-a2a4-9372bea9f5aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import glob\n",
    "import sys\n",
    "import shutil\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "from mintpy.utils import readfile, writefile, utils as ut\n",
    "from mintpy import workflow, view, tsview, plot_network, plot_transection, plot_coherence_matrix\n",
    "\n",
    "from mintpy.objects import timeseries\n",
    "from datetime import datetime\n",
    "import numpy as np\n",
    "from skimage.morphology import remove_small_objects\n",
    "\n",
    "from scipy.stats import linregress\n",
    "import pandas as pd\n",
    "import matplotlib\n",
    "import datetime"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "694f1ed6-f6b9-4ca6-9679-c59183b22d9c",
   "metadata": {},
   "source": [
    "['load_data', 'modify_network', 'reference_point', 'quick_overview', 'correct_unwrap_error']\n",
    "  ['invert_network', 'correct_LOD', 'correct_SET', 'correct_troposphere', 'deramp', 'correct_topography']\n",
    "  ['residual_RMS', 'reference_date', 'velocity', 'geocode', 'google_earth', 'hdfeos5']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1902a919-9d84-4c41-ba8c-24c8f8de0ad1",
   "metadata": {},
   "outputs": [],
   "source": [
    "cd /home/micromamba/data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d81d5ead-f4be-4da8-b77e-af16dc9e9de3",
   "metadata": {},
   "outputs": [],
   "source": [
    "cp /home/micromamba/scripts/configs/desc_upper.cfg ./hyp3.cfg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd5dd34b-675c-482a-b2ae-b1db9191b924",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## load_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ddc86f7-f20e-4b74-8c07-112db3f279b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "!load_data.py -t hyp3.cfg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "000bd76c-b9f6-4845-9c03-3eef38082abb",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## double_difference\n",
    "\n",
    "This the only step that isn't included in the normal MintPy `smallbaselineApp` workflow. It highlights local deformation, but obscures the presence of large-scale (greater than 10 km scale) deformation and unwrapping errors. See [Bekaert et. al 2020](https://www.sciencedirect.com/science/article/pii/S0034425720303539) for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faf27801-c5ad-418d-ad26-efa65170c48d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "desired_width = 2000\n",
    "pixel_size = 40\n",
    "\n",
    "x = 750\n",
    "n_pixels = np.floor((desired_width / 40) /2)\n",
    "width = ((n_pixels*2)+1)*400\n",
    "print(f'N pixels: {n_pixels}\\nActual width: {width}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5c0e360-0e63-4ac3-9525-e408926db694",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from skimage import morphology\n",
    "filter_par = [2,25]\n",
    "\n",
    "local_kernel = morphology.disk(filter_par[0], np.float32)\n",
    "regional_kernel = morphology.disk(filter_par[1], np.float32)\n",
    "\n",
    "center = slice(filter_par[1]-filter_par[0],filter_par[1]+filter_par[0]+1)\n",
    "print(center,center)\n",
    "regional_kernel[center,center] =  regional_kernel[center,center] - local_kernel\n",
    "\n",
    "local_kernel = local_kernel / local_kernel.flatten().sum()\n",
    "# local_filt = ndimage.convolve(data, local_kernel)\n",
    "\n",
    "regional_kernel = regional_kernel / regional_kernel.flatten().sum()\n",
    "# regional_filt = ndimage.convolve(data, regional_kernel)\n",
    "\n",
    "f, (ax1,ax2) = plt.subplots(1,2)\n",
    "img1 = ax1.imshow(local_kernel,interpolation='none')\n",
    "ax1.set_title(f'Width: {local_kernel.shape[0]*40}m')\n",
    "img2 = ax2.imshow(regional_kernel,interpolation='none')\n",
    "ax2.set_title(f'Width: {regional_kernel.shape[0]*40}m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f84637fd-e3fb-4fec-b06e-6ae3c72632b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "params = [2,25]\n",
    "outname = f'inputs/ifgramStack_{params[0]}_{params[1]}.h5'\n",
    "!spatial_filter.py inputs/ifgramStack.h5 unwrapPhase -o {outname} -f double_difference -p {params[0]} {params[1]}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4123fac7-8a46-4016-879d-ac9ae09e8cde",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## directory management"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f04c7ca8-ee2d-48a9-86a3-63154396e3de",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6644dab9-a76f-4b7c-9b46-073ecb03e715",
   "metadata": {},
   "outputs": [],
   "source": [
    "mkdir inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cacde1e-683d-46eb-9b6d-64f917bc97b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "cp -r /home/micromamba/raw/desc_upper/mintpy/inputs/ifgramStack_2_25.h5  ./inputs/ifgramStack_2_25.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "388df030-ff1d-42f4-abdd-9d0f29e54237",
   "metadata": {},
   "outputs": [],
   "source": [
    "cp -r /home/micromamba/raw/desc_upper/mintpy/inputs/geometryGeo.h5  ./inputs/geometryGeo.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c563eff4-97d2-45eb-8505-ec684bcd7781",
   "metadata": {},
   "outputs": [],
   "source": [
    "# rm -rf inputs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac271f65-001e-4feb-b7ca-90419cc261ee",
   "metadata": {
    "tags": []
   },
   "source": [
    "## generate water mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a0921b2-fb55-495c-b207-dfb72b949e6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "!generate_mask.py ./inputs/geometryGeo.h5 waterMask -m 0.5 -o waterMask.h5\n",
    "view.main('waterMask.h5 --noverbose'.split())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72d94f95-2f7e-4d58-91c6-968cf3c08945",
   "metadata": {
    "tags": []
   },
   "source": [
    "## reference_point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c10f9b1-1a2d-47ec-a7e1-db125f814259",
   "metadata": {},
   "outputs": [],
   "source": [
    "view.main('inputs/geometryGeo.h5 height -m waterMask.h5 --noverbose'.split())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38706a0c-17b4-48d2-93a4-9f13c4a1a5d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "!temporal_average.py inputs/ifgramStack_2_25.h5 -d coherence -o avgSpatialCoh.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07417891-7f07-4cd8-804c-95cfa5019540",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "view.main('avgSpatialCoh.h5 --noverbose --figsize 8 8'.split())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a1c6351-97a9-44c2-90b5-f24f9f608d64",
   "metadata": {},
   "outputs": [],
   "source": [
    "!reference_point.py inputs/ifgramStack_2_25.h5 --reset\n",
    "!reference_point.py inputs/ifgramStack_2_25.h5 -t ./hyp3.cfg -m waterMask.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a89812b4-8abd-4298-8499-1dc586adc30f",
   "metadata": {
    "tags": []
   },
   "source": [
    "## modify_network"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac21a1b6-6039-4106-9b69-7b042415ee6e",
   "metadata": {},
   "source": [
    "This step has the largest effect on the overall quality of the analysis so **take your time here and remove garbage data!**\n",
    "Remember the wise words:\n",
    "> Garbage in, Garbage out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a94b3e4-8659-4cff-b72d-c08d30f8ae83",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "!spatial_average.py inputs/ifgramStack.h5 -d coherence -m waterMask.h5\n",
    "\n",
    "m = pd.read_csv('coherenceSpatialAvg.txt',skiprows=5,header=None,sep='\\t',\n",
    "                names=['DATE12','Mean','Btemp/days','Bperp/m','Num'])\n",
    "\n",
    "m['date1'],m['date2'] = zip(*m['DATE12'].str.split('_'))\n",
    "m['date1'] = pd.to_datetime(m['date1'])\n",
    "m['date2'] = pd.to_datetime(m['date2'])\n",
    "\n",
    "dates = sorted(list(set(list(m.date1) + list(m.date1))))\n",
    "\n",
    "matrix = pd.DataFrame(index=dates,columns=dates)\n",
    "for i,row in m.iterrows():\n",
    "    matrix.loc[row['date1'],row['date2']] = row['Mean']\n",
    "    matrix.loc[row['date2'],row['date1']] = row['Mean']\n",
    "    \n",
    "avg_coh = pd.DataFrame({'date':matrix.index,'coherence':matrix.mean(axis=1)}).reset_index(drop=True)\n",
    "\n",
    "smooth_coh = avg_coh.copy()\n",
    "smooth_coh.index = pd.DatetimeIndex(smooth_coh.date)\n",
    "smooth_coh = smooth_coh.sort_index()\n",
    "smooth_coh = smooth_coh.rolling(datetime.timedelta(days=60),min_periods=1).mean()\n",
    "smooth_coh = smooth_coh.reset_index()\n",
    "\n",
    "f, ax = plt.subplots(1,1,figsize=(15,5))\n",
    "\n",
    "ax.plot(avg_coh.date,avg_coh.coherence)\n",
    "ax.plot(smooth_coh.date,smooth_coh.coherence,color='red')\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Coherence')\n",
    "plt.savefig('coh_timeseries.png')\n",
    "# plt.xlim((datetime.datetime(2021,3,1),datetime.datetime(2021,4,1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de26bf6b-ad26-4f64-b614-46720162ed00",
   "metadata": {},
   "source": [
    "Based on this graph, there is a clear seasonal trend in the correlation data. I've found that performing the SBAS inversion on each year separately (breaking the analysis at the low points) produces the best results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5289472f-c207-4b1c-baf4-e6b37b28c280",
   "metadata": {},
   "outputs": [],
   "source": [
    "cp /home/micromamba/scripts/configs/desc_upper.cfg ./hyp3.cfg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "518c8f6e-6747-488c-bfeb-cd0e8a4eb6d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "!modify_network.py inputs/ifgramStack_2_25.h5 --reset\n",
    "!modify_network.py inputs/ifgramStack_2_25.h5 -t hyp3.cfg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13e74d80-58d1-4d60-9d75-c2ebdf8326a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_network.main('inputs/ifgramStack_2_25.h5 -t hyp3.cfg --show-kept --figsize 15 5'.split())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92c789ec-c4c2-4aea-bba4-1ecf806c7a65",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## view ifgram stack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e82b4a96-caba-458b-bedc-8e375c0777e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = readfile.read_template('hyp3.cfg')\n",
    "start_date = datetime.datetime.strptime(d['mintpy.network.startDate'], '%Y%m%d')\n",
    "end_date = datetime.datetime.strptime(d['mintpy.network.endDate'], '%Y%m%d')\n",
    "refy, refx = d['mintpy.reference.yx'].split(', ')\n",
    "min_coh = float(d['mintpy.network.minCoherence'])\n",
    "\n",
    "used_num = [str(x) for x in m.loc[(m.date1 >= start_date) & (m.date2 <= end_date), 'Num']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28c52f6b-b6aa-46a4-9187-7b16e8ac5eb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "view.main(f'inputs/ifgramStack_2_25.h5 --noverbose --ref-yx {refy} {refx} -n {\" \".join(used_num)} --vlim -3 3'.split())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b25389ca-e0d3-41d8-8015-9bde4ed2660d",
   "metadata": {
    "tags": []
   },
   "source": [
    "## invert_network\n",
    "This can take a while so write it to a log if you need to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52fdb4ea-ea15-4b8a-aa96-e74c1cd37809",
   "metadata": {},
   "outputs": [],
   "source": [
    "# !ifgram_inversion.py inputs/ifgramStack.h5 -t ./hyp3.cfg |& tee -i invert.log ; python ../JPL_internship/emailme.py\n",
    "!ifgram_inversion.py inputs/ifgramStack_2_25.h5 -t ./hyp3.cfg -m waterMask.h5 -o timeseries_descupper_2021.h5 temporalCoherence_descupper_2021.h5 numInvIfgram_descupper_2021.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2205a97c-ce59-4334-a576-bfd88daf8474",
   "metadata": {
    "tags": []
   },
   "source": [
    "## check temporal ramp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4377e7c8-cf61-41fa-b673-d84ae97a90c9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ts = timeseries(file='timeseries_descupper_2021.h5')\n",
    "\n",
    "data = ts.read()\n",
    "data[data == 0] = np.nan\n",
    "# data[:,~coh_mask] = np.nan\n",
    "dates = pd.to_datetime(ts.get_date_list())\n",
    "means = np.nanmean(data,axis=(1, 2))\n",
    "means[0] = 0\n",
    "# datetimes = np.array([datetime.strptime(x,'%Y%m%d').toordinal() for x in dates])\n",
    "\n",
    "datetimes = np.array([int(x.toordinal()) for x in dates])\n",
    "slope, intercept, r, p, se = linregress(datetimes,means)\n",
    "fit = slope*datetimes+intercept\n",
    "print(f'slope: {round(slope*100,4)} cm\\nr2:{round(r**2,2)}')\n",
    "\n",
    "plt.scatter(datetimes,means)\n",
    "plt.plot(datetimes,fit,color='red')\n",
    "plt.ylim(-0.02,0.02)\n",
    "plt.xlabel('Ordinal Date')\n",
    "plt.ylabel('Avg Timeseries Value')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c518500-7e85-4024-913d-f437d948fe22",
   "metadata": {},
   "source": [
    "## velocity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b0ba9a3-6b81-4f8e-9271-613aaa0daeba",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "outname = 'velocity_descupper_2021.h5'\n",
    "!timeseries2velocity.py timeseries_descupper_2021.h5 -o {outname} --bootstrap --bc 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72a6bbd3-cb32-4373-b38d-a2d964cf169b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "view.main(f'{outname} --noverbose -m no --vlim -5 5'.split())\n",
    "# view.main('temporalCoherence_descupper_2021.h5 --noverbose --colormap RdBu_truncate --cm-vlist 0.0 0.7 1.0'.split())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe920cbf-ed21-4fbc-b917-573dca8065b7",
   "metadata": {},
   "source": [
    "## create mask file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe85b2e5-7f3d-467c-bf25-fbbe84579638",
   "metadata": {},
   "outputs": [],
   "source": [
    "!generate_mask.py temporalCoherence_descupper_2021.h5 -m 0.01 -o maskTempCoh_descupper.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "765aada7-47b9-4399-bc43-bbf132035355",
   "metadata": {},
   "outputs": [],
   "source": [
    "view.main(f'maskTempCoh_descupper.h5 --noverbose'.split())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2dd8e65-ea36-4086-b450-cb2cdd18c3ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "ls *h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9bca91a-8c4d-4668-b801-26178a93b146",
   "metadata": {},
   "outputs": [],
   "source": [
    "cp *h5 /mnt/i/hyp3_downloads/desc_upper/mintpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "176f23a6-7f2b-4233-9df8-4552df345389",
   "metadata": {},
   "source": [
    "## complete!"
   ]
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}