full_data_downsample_kite.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Downsampling data with <i>kite</i>\n",
    "\n",
    "The <i><a href=\"https://pyocko.org\">pyrocko</a></i> suite of tools comprise a series of open source applications designed to simplify the modeling of earthquakes, using both seismic and InSAR data. Here we will use one of those applications, <i>kite</i>, to read in and downsample an ISCE interferogram using quadtree decomposition.\n",
    "\n",
    "## Setup (optional)\n",
    "\n",
    "First, install some dependencies, if you haven't already..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the kite and pyrocko codes - ignore if you already have these installed!\n",
    "!pip install pyrocko\n",
    "!pip install utm\n",
    "!pip install git+https://github.com/pyrocko/kite"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...and import the necessary dependencies: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from kite import Scene\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm, colors\n",
    "import math\n",
    "from osgeo import gdal\n",
    "import utm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we specify a path and some filenames:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ifgpath=\"/home/gareth/work/quakes/elazig/merged/\"                     # directory with your interferogram\n",
    "ifgname=\"filt_topophase.unw.geo\"                                      # interferogram file\n",
    "corname=\"topophase.cor.geo\"                                           # correlation file\n",
    "lospath=ifgpath                                                       # can specify different path for los information\n",
    "losname=\"los.rdr.geo\"                                                 # line-of-sight file\n",
    "mskname=\"water.msk\"                                                   # water mask file (optional)\n",
    "use_water_mask=True\n",
    "\n",
    "outbase=\"elazig_asc\"                                                  # output filename base"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are a couple of important parameters. <i>kite</i> currently cannot read in the radar wavelength from the interferogram metadata, as it is not included in the xml file(!) and so we specify it here. We also supply a correlation threshold for masking out low correlation data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wavel = 0.0555     # radar wavelength in m\n",
    "cor_thresh = 0.25  # correlation threshold between 0 (bad) and 1 (perfect)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load in the interferogram\n",
    "\n",
    "<i>kite</i> has a single command to read in an interferogram from ISCE (see below), but manual import is more flexible and powerful, allowing for cropping and masking. We make use of all of that functionality here. The steps are:\n",
    "\n",
    "- Read in the interferogram, correlation and LOS files (and optionally apply a water mask)\n",
    "- Plot the input interferogram in image coordinates\n",
    "- Specify crop bounds in image coordinates if you want\n",
    "- Convert to displacement in meters and apply a correlation mask\n",
    "- Extract coordinate system metadata from the imported images\n",
    "- Crop the images and import them into a <i>kite</i> scene object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# that single command, for reference\n",
    "#sc = Scene.import_data(ifgpath+ifgname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming we want to do everything manually, first we read in the interferogram using gdal, and plot it up with matplotlib to see if we want to crop it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load in your data manually!\n",
    "\n",
    "# unwrapped interferogram\n",
    "ifgfile=gdal.Open(ifgpath+ifgname,gdal.GA_ReadOnly)   # open with gdal\n",
    "rb = ifgfile.GetRasterBand(2)                         # phase info is in band 2\n",
    "ifg = rb.ReadAsArray()                                # read it in as a number array\n",
    "\n",
    "# let's plot it to have a look\n",
    "### matplotlib inline\n",
    "fig, ax = plt.subplots(figsize=(10,10))        # new figure called \"ax1\"\n",
    "im = ax.imshow(ifg,origin='upper')             # plot displacements with origin at upper left (ISCE default)\n",
    "fig.colorbar(im)                               # plot a color bar!\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to crop it, specify the bounds here. And we will apply it, and all of the other steps, in the code below! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# specify crop bounds (in image coordinates, e.g. from mdx, in pixels from the top left)\n",
    "# for the full scene, set minima to 0 and maxima to -1\n",
    "ymin=0\n",
    "ymax=-1\n",
    "xmin=0\n",
    "xmax=-1\n",
    "\n",
    "ifg=np.flip(ifg,0)                                    # flip your interferogram to account for origin position\n",
    "                                                      # (kite prefers lower left; ISCE, top left)\n",
    "# correlation file\n",
    "corfile=gdal.Open(ifgpath+corname,gdal.GA_ReadOnly)   # open with gdal\n",
    "rb = corfile.GetRasterBand(2)                         # correlation info is in band 2\n",
    "cor = rb.ReadAsArray()                                # read it in as a number array\n",
    "cor=np.flip(cor,0)                                    # flip image to account for origin position\n",
    "                                                      # (kite prefers lower left; ISCE, top left)\n",
    "# line-of-sight information\n",
    "losfile=gdal.Open(lospath+losname,gdal.GA_ReadOnly)   # open with gdal\n",
    "rb = losfile.GetRasterBand(1)                         # incidence info is in band 1\n",
    "inc = rb.ReadAsArray()                                # read it in as a number array\n",
    "inc=np.flip(inc,0)\n",
    "rb = losfile.GetRasterBand(2)                         # phase info is in band 2\n",
    "azi = rb.ReadAsArray()                                # read it in as a number array\n",
    "azi=np.flip(azi,0)   \n",
    "    \n",
    "# apply the water mask, if you have one\n",
    "if use_water_mask==True:\n",
    "    # water mask file\n",
    "    mskfile=gdal.Open(ifgpath+mskname,gdal.GA_ReadOnly)  \n",
    "    rb = mskfile.GetRasterBand(1)                         # only one band here\n",
    "    msk = rb.ReadAsArray()\n",
    "    msk=np.flip(msk,0)    \n",
    "    ifg[msk<1]=np.nan                                     # mask out water\n",
    "\n",
    "\n",
    "# convert displacements to meters and mask low correlation pixels\n",
    "ifg*=wavel/4/math.pi                # change values of sc.displacement in place\n",
    "ifg[cor<cor_thresh]=np.nan          # mask out low correlation\n",
    "\n",
    "# get metadata for the image: frame size, pixel size and reference point\n",
    "gt = losfile.GetGeoTransform()    # geotransform information\n",
    "nrows=losfile.RasterYSize\n",
    "\n",
    "# sort out the maximum crop size in the y direction\n",
    "if ymax==-1:\n",
    "    ymax=nrows\n",
    "\n",
    "# now for the kite part...\n",
    "    \n",
    "# initiate your kite scene object\n",
    "sc = Scene()\n",
    "    \n",
    "# specify pixel spacing information\n",
    "sc.frame.spacing = 'degree'\n",
    "sc.frame.dN = abs(gt[5])\n",
    "sc.frame.dE = gt[1]\n",
    "\n",
    "# set your reference point\n",
    "sc.frame.llLat = gt[3]+ymax*gt[5]\n",
    "sc.frame.llLon = gt[0]+xmin*gt[1]\n",
    "\n",
    "# read your interferogram and LOS information into the scene object (and crop!)\n",
    "sc.displacement=ifg[nrows-ymax:nrows-ymin,xmin:xmax]\n",
    "sc.theta=np.radians(90-inc[nrows-ymax:nrows-ymin,xmin:xmax])\n",
    "sc.phi=np.radians(90+azi[nrows-ymax:nrows-ymin,xmin:xmax])\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can look and see what we imported into <i>kite</i> using matplotlib..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### matplotlib inline\n",
    "fig, ax = plt.subplots(figsize=(10,10))        # new figure called \"ax1\"\n",
    "im = ax.imshow(sc.displacement,origin='lower') # plot displacements with origin at lower left (kite default)\n",
    "fig.colorbar(im)                                # plot a color bar!\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quadtree decomposition\n",
    "\n",
    "We can initialize the quadtree decomposition in <i>kite</i> by defining an abbreviation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qt = sc.quadtree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can specify some parameters that control the quadtree process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qt.epsilon = 0.01     # Variance threshold\n",
    "qt.nan_allowed = 0.1      # Percentage of NaN values allowed per tile/leaf\n",
    "qt.tile_size_max = 0.1    # Maximum leaf size in degrees\n",
    "qt.tile_size_min = 0.01  # Minimum leaf size in degrees"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the distribution of leaves/tiles. Typically we are looking for a few hundred... If you don't like what you have, try editing the values above and re-running. If you would like more leaves, lower the epsilon value and/or make the tile sizes smaller. If you want fewer, do the opposite!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"There are\",qt.nleaves,\"leaves\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And if you want to see the quadtree results, you can. This example prints out the leaf number, the lower left ('ll') east and north coordinates of the leaf (with respect to the lower left of the whole image), and the mean displacement of the pixels in that leaf:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "j=0\n",
    "for k in qt.leaves:\n",
    "    print(\"leaf\",j,\"llE:\",qt.leaves[j].llE,\"llN:\",qt.leaves[j].llN,\"mean:\",qt.leaves[j].mean) \n",
    "    j=j+1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And you can plot them too!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matplotlib inline\n",
    "fig, ax3 = plt.subplots(figsize=[5,7.5])\n",
    "\n",
    "ulimit = qt.leaf_medians.max()\n",
    "llimit = qt.leaf_medians.min()\n",
    "\n",
    "color_map = cm.ScalarMappable(\n",
    "    norm=colors.Normalize(vmin=llimit, vmax=ulimit))\n",
    "\n",
    "for rect, leaf in zip(qt.getMPLRectangles(), qt.leaves):\n",
    "    color = color_map.to_rgba(leaf.median)\n",
    "    rect.set_facecolor(color)\n",
    "    ax3.add_artist(rect)\n",
    "\n",
    "\n",
    "ax3.set_xlim(0, qt.leaf_eastings.max())\n",
    "ax3.set_ylim(0, qt.leaf_northings.max())\n",
    "ax3.set_aspect('equal')\n",
    "\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Processing the quadtree output\n",
    "\n",
    "There are several modeling codes out there for InSAR data, and there is no standard data input format. But typically, all such modeling codes will need (at a minimum) data point locations, displacements and line-of-sight vector information.\n",
    "\n",
    "Here is an example of reformatting the output of the <i>kite</i> quadtree decomposition for use in one such code, <i>okinv</i>. This involves extracting the quadtree leaf centroids and displacement averages in their lat-long coordinates, and then converting them to a local coordinate system (a UTM system, with coordinates specified in kilometers). We also include the LOS vector information for each of these points. Hopefully this is straightforward enough to follow that you can reformat the data for use in other codes, if necessary!\n",
    "\n",
    "The steps involved here:\n",
    "\n",
    "- Convert the reference point from lat-long coordinates to UTM\n",
    "- Compute the three components of the unit LOS vector for each pixel\n",
    "- Convert the positions of the rest of the points to UTM, and plot them\n",
    "\n",
    "The reference point in <i>kite</i> is the lower left (SW) corner of the interferogram. (Note that <i>ISCE</i>, and indeed most geospatial imaging formats, use the upper left corner.) We can use the <i>utm</i> code to convert them to rectilinear coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reflat=sc.frame.llLat\n",
    "reflon=sc.frame.llLon\n",
    "\n",
    "# we would also like the reference point in utm coordinates\n",
    "refutm=utm.from_latlon(reflat,reflon)\n",
    "\n",
    "print(\"reference point (lat-long):\",reflat,\"N\",reflon,\"E\")\n",
    "\n",
    "print(\"reference point (utm):\",refutm[0],\"E\",refutm[1],\"N\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extracting the information we need from <i>kite</i>'s output is fairly straightforward. 'leaf_focal_points' is the name given to the centers of each quadtree cell, and is provided both in meters in a local coordinate system, as well as in long-lat.\n",
    "\n",
    "The definitions of incidence ('qt.leaf_thetas') and azimuth ('qt.leaf_phis') in kite are a little unusual (as seems to be common in InSAR, everyone has their own convention!) Here, the convention is up-looking, towards the satellite, with the azimuth being counter-clockwise from the East. Both values are in radians.\n",
    "\n",
    "We can convert this fairly easily to other definitions such as 'range change', i.e. with the unit line-of-sight (LOS) vector pointing from the satellite to the ground target, by doing some simple trigonometry. (It is probably a good idea to draw the problem, to see how this works...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# xy positions of data points in long-lat\n",
    "xyll=np.column_stack((qt.leaf_focal_points[:,0]+reflon,qt.leaf_focal_points[:,1]+reflat))\n",
    "\n",
    "# los vector components using kite definiton of coordinates\n",
    "losx=np.cos(qt.leaf_thetas)*-np.cos(qt.leaf_phis)\n",
    "losy=np.cos(qt.leaf_thetas)*-np.sin(qt.leaf_phis)\n",
    "losz=-np.sin(qt.leaf_thetas)\n",
    "\n",
    "# okinv format includes an id label for each datapoint\n",
    "dataptid=np.arange(1,losx.size+1)\n",
    "\n",
    "# and finally let's convert our lat-long coordinates to UTM km\n",
    "xym = utm.from_latlon(xyll[:,1],xyll[:,0]) \n",
    "xkm = xym[0]/1000\n",
    "ykm = xym[1]/1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the points, just to make sure that our conversion makes sense:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,10))\n",
    "\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_xlabel(\"UTM_x (km)\",fontsize=12)\n",
    "ax.set_ylabel(\"UTM_y (km)\",fontsize=12)\n",
    "ax.grid(True,linestyle='-',color='0.75')\n",
    "\n",
    "# scatter with colormap mapping to z value\n",
    "im=ax.scatter(xkm,ykm,s=20,c=qt.leaf_means, marker = 'o');\n",
    "#im=ax.scatter(xkm,ykm,s=20,c=qt.leaf_means, marker = 'o', cmap = cm.jet); # if you prefer jet\n",
    "ax.set_aspect('equal')\n",
    "fig.colorbar(im) \n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving the output\n",
    "\n",
    "We finish by outputting the data to file, in one example format. The following format, based on the okinv code, has seven columns:\n",
    "\n",
    "\\<x_coordinate(km)\\>  \\<y_coordinate(km)\\>  \\<los_displacement(m)\\>  \\<unit_los_x\\>  \\<unit_los_y\\>  \\<unit_los_z\\>  \\<id\\>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# okinv format in km\n",
    "okinv_data=np.column_stack((xkm,ykm,qt.leaf_means,losx,losy,losz,dataptid))\n",
    "outfile=outbase+\".okinv\"\n",
    "np.savetxt(outfile,okinv_data,fmt='%f %f %f %f %f %f %d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, if your code prefers it (and ours does not), you can save the data in long-lat coordinates:\n",
    "\n",
    "\\<longitude\\> \\<latitude\\> \\<los_displacement(m)\\> \\<unit_los_x\\> \\<unit_los_y\\> \\<unit_los_z\\> \\<id\\>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# okinv-equivalent format in lat-long coordinates \n",
    "okinv_data_ll=np.column_stack((xyll,qt.leaf_means,losx,losy,losz,dataptid))\n",
    "outfile=outbase+\"_ll\"+\".okinv\"\n",
    "np.savetxt(outfile,okinv_data_ll,fmt='%f %f %f %f %f %f %d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also, optionally, save the data in <i>kite</i>'s internal format, and (potentially, per the <i>kite</i> manual) in a comma-separated variable file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save the scene information in a single file\n",
    "sc.save(outbase)\n",
    "# export the quadtree information in a CSV file\n",
    "# qt.export(outbase+'.csv')   # not working for me right now..."
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}