test merging overlapping images with gdal python pixel function

pixel-function-test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test gdal vrt pixel functions\n",
    "\n",
    "Use GDAL python pixel function to define how overlapping image regions are dealt with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import rasterio\n",
    "from rasterio.transform import Affine\n",
    "import rasterio.plot\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create two simple rasters, create vrts overlapping by 1 degree, use pixel functions to mosaiic\n",
    "x = np.linspace(-4.0, 4.0, 240)\n",
    "y = np.linspace(2.0, 8.0, 180)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "Z1 = np.ones_like(X)\n",
    "\n",
    "Z1[:10,:] = 0 #or np.nan?\n",
    "Z1[:,:10] = 0\n",
    "\n",
    "Z1c = Z1.astype('complex64') + 1j\n",
    "\n",
    "res = (x[-1] - x[0]) / 240.0\n",
    "transform1 = Affine.translation(x[0] - res / 2, y[-1] - res / 2) * Affine.scale(res, -res)\n",
    "transform1\n",
    "\n",
    "with rasterio.open('top-complex.tif', 'w',\n",
    "                    driver='GTiff',\n",
    "                    height=Z1c.shape[0],\n",
    "                    width=Z1c.shape[1],\n",
    "                    #nodata=np.nan, #0 #throwing KeyError: 'complex64'... likely a bug\n",
    "                    count=1,\n",
    "                    dtype=Z1c.dtype,\n",
    "                    crs='+proj=latlong',\n",
    "                    transform=transform1,\n",
    "                    # creation options\n",
    "                    blockxsize=128,\n",
    "                    blockysize=128,\n",
    "                    tiled=True,\n",
    "                    ) as dst:\n",
    "    dst.write(Z1c, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create two simple rasters, create vrts overlapping by 1 degree, use pixel functions to mosaiic\n",
    "x = np.linspace(-4.0, 4.0, 240)\n",
    "y = np.linspace(-3.0, 3.0, 180)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "Z2 = 2*np.ones_like(X)\n",
    "\n",
    "Z2[-10:,:] = 0\n",
    "Z2[:,-10:] = 0\n",
    "\n",
    "Z2c = Z2.astype('complex64') + 2j\n",
    "\n",
    "res = (x[-1] - x[0]) / 240.0\n",
    "transform2 = Affine.translation(x[0] - res / 2, y[-1] - res / 2) * Affine.scale(res, -res)\n",
    "transform2\n",
    "\n",
    "with rasterio.open('bottom-complex.tif', 'w',\n",
    "                    driver='GTiff',\n",
    "                    height=Z2c.shape[0],\n",
    "                    width=Z2c.shape[1],\n",
    "                    #nodata=np.nan,\n",
    "                    count=1,\n",
    "                    dtype=Z2c.dtype,\n",
    "                    crs='+proj=latlong',\n",
    "                    transform=transform2,\n",
    "                    # creation options\n",
    "                    blockxsize=128,\n",
    "                    blockysize=128,\n",
    "                    tiled=True,\n",
    "                    ) as dst:\n",
    "    dst.write(Z2c, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'driver': 'GTiff', 'dtype': 'complex64', 'nodata': None, 'width': 240, 'height': 180, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.03333333333333333, 0.0, -4.016666666666667,\n",
      "       0.0, -0.03333333333333333, 7.983333333333333), 'blockxsize': 128, 'blockysize': 128, 'tiled': True, 'interleave': 'band'}\n",
      "{'driver': 'GTiff', 'dtype': 'complex64', 'nodata': None, 'width': 240, 'height': 180, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.03333333333333333, 0.0, -4.016666666666667,\n",
      "       0.0, -0.03333333333333333, 2.9833333333333334), 'blockxsize': 128, 'blockysize': 128, 'tiled': True, 'interleave': 'band'}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-3, 8)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD8CAYAAACxd9IeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAWGklEQVR4nO3de5AdZZnH8e+Pc2ZA7kgghCQo61IWeGHFFODGWsALBgQtqywXVHBRK6ulrlpe1kutt/3HKnZdtRDdFLIsBWi5SpTaighrueVtoQwRLxiwsiELuZNECAg6c2ae/aM78TCZOaffOZ3pPtO/T9WpmdP99ttvUvXM2/322++jiMDM5r9Dqm6Amc0NB7tZQzjYzRrCwW7WEA52s4ZwsJs1RCnBLun9ku6T9GtJX5N0WBn1ms1XkpZK+oGk9XnsvHeaMpL0RUkbJP1S0lld+1ZIeiDf95Ei5xw42CUtBv4OWBYRzwdawGWD1ms2z3WAD0TE6cC5wLsknTGlzEXAaflnJfBlAEkt4Ev5/jOAy6c59gBlXca3gWdIagOHA1tLqtdsXoqIbRGxLv/9cWA9sHhKsdcCN0bmLuBYSYuAs4ENEbExIsaAr+dle2qX0Ogtkv4JeAh4CrgjIu6YWk7SSrK/Thw6cuiLFy2Y+u8yO7g2bdu4KyJOGKSOV73i9Ni95/eFyt7z84fvA/7QtWlVRKyaWk7Ss4EXAXdP2bUYeLjr++Z823Tbz+nXnoGDXdJxZH9VTgUeBf5D0psj4qbucvk/chXAqSc/Jz7zt58d9NRmSa781Bv+b9A6du/ayd3fX1GobHvBLX+IiGW9ykg6EvgW8L6I2Dt19zSHRI/tvdvTr0ABrwAejIhHACTdCvwlcFPPo8yGUMQYE52HSqlL0ghZoN8cEbdOU2QzsLTr+xKyW+TRGbb3VEawPwScK+lwssv4lwNrS6jXrH40wmTr5MGrkQR8FVgfEZ+bodhtwLslfZ3sMv2xiNgm6RHgNEmnAlvIBsTf2O+cZdyz3y3pm8A6shHGn5NfrpvNNwFM9r9iLmI5cAXwK0n35ts+BpwCEBFfAdYAFwMbgCeBq/J9HUnvBr5H9vTr+oi4r98Jy+jZiYhPAp8soy6zupss4bXwiPgx0997d5cJ4F0z7FtD9segsFKC3awpgiirZ59ztQr2Pa/y4zg7iD5VTjWT5VQz52oV7GZ1FzHG+EQ5o/FzzcFulkKj0FpSsPDPDmpTUjnYzRIEMDGct+wOdrNUvmc3a4gh7dgd7GYpJmOMP3Q2V92MWXGwmyUZ5ZDW0v7FAA/QmQ25yd4T32rLwW6WREyGg91s3stehHGwmzWCg92sASZjjCfHt1TdjFlxsJslkA6l3T6lYOmfHNS2pHKwmyXI7tmHM7eKg90ske/ZzRogwI/ezJpgMsZ5olNODhRJ1wOXADvzbEpT938IeFP+tQ2cDpwQEXskbQIeByaATr8lq/dVYGYFSaOMtooO0PV1A3ANcON0OyPiauDq7Ly6FHh/ROzpKnJBROwqejIHu1misu7ZI+KHeTaYIi4HvjbI+YZzWNGsIoGYLPgpS56TYQVZQok/NQXukHRPnlqtL/fsZikCJqNwH7lAUnfClGlzvRVwKfCTKZfwyyNiq6QTgTsl3R8RP+xViYPdLFFCr72ryMBZAZcx5RI+IrbmP3dKWk2W2dXBblaWiRhn7/j2OTufpGOA84A3d207AjgkIh7Pf78Q+Ey/uhzsZgmkUQ5rF11dtl9d+hpwPtnl/mayrEojsD/9E8DryNKgd+eJXgisztLF0QZuiYjb+53PwW6WqKxJNRFxeYEyN5A9ouvethE4M/V8DnazRJ4ua9YA2aO34Xxi7WA3S+S58WYN0IlxHh3fUXUzZqWUYJd0LHAd8HyymT1vjYj/KaNuszppaZTD28OZbbisnv0LwO0R8XpJo8DhJdVrViuNXnBS0tHAXwF/AxARY8DYoPWa1VJADOk9exnDin8GPAL8m6SfS7oun9XzNJJWSlorae3jT+4t4bRm1ZjrF2HKUsZlfBs4C3hPRNwt6QvAR4B/6C6UvwCwCuDUk58zrLnxrOE60WHP2CNVN2NWygj2zcDmiLg7//5NsmA3m3cO0QhHjiyquhmzMvBlfERsBx6W9Nx808uB3wxar1ldTYYKfeqmrNH49wA35yPxG4GrSqrXrFai6bneIuJeoIz3ds1qr46Db0V4Bp1ZonCwm81/nckOu/5YeEHXWnGwmyVoaYRjRhZW3YxZcbCbJXBGGLMG8QCdWUMMa88+nEtumFUkyEbji3z6kXS9pJ2Sfj3D/vMlPSbp3vzzia59KyQ9IGmDpEIzVt2zmyXoxAQ7/7i7rOpuoEeut9yPIuKS7g2SWsCXgFeSTVf/maTbIqLnzFUHu1mCltocN3JiKXUl5nrrdjawIV9lFklfB15Ln2nqvow3S5TwiuuCfa91559COdmmeImkX0j6rqTn5dsWAw93ldmcb+vJPbtZirTFKwZN/7QOeFZEPCHpYuDbwGkw7YBA39fG3bObJZjLLK4RsTcinsh/XwOMSFpA1pMv7Sq6BNjarz737GYJOtFhx1N7+hcsgaSTgB0REZLOJuucdwOPAqdJOhXYQpb48Y396nOwmyVoq83xh55QSl0Fcr29HninpA7wFHBZRATQkfRu4HtAC7g+Iu7r2/ZSWm3WENl02ZLq6pPrLSKuIXs0N92+NcCalPM52M0S+RVXs0Zo+Eo1Zk0RQ7xuvIPdLMF4dNj2h7kZjS+bg90sQVttFhy6oOpmzIqD3SyRL+PNGqKsR29zzcFulqTYu+p15GA3SxD4Mt6sEcYnO2x96ndVN2NWHOxmCUbU5sRDj6+6GbPiYDdL4Mt4swYZ0sF4B7tZGrlnN2uCsckOW558tOpmzEppwZ4vb7sW2DJ16Vuz+WJEbU467JlVN2NWylyD7r3A+hLrM6ulKPipm1KCXdIS4NXAdWXUZ1ZnESr0qZuyevbPAx8GJmcqIGnlvvWzH39yb0mnNZt72Tvt/T/9FEj/9CZJv8w/P5V0Zte+TZJ+laeFWluk3QMHu6RLgJ0RcU+vchGxKiKWRcSyow4/etDTmlWmxJ79BmBFj/0PAudFxAuBfwRWTdl/QUT8RdG16csYoFsOvCZfxP4w4GhJN0XEm0uo26xWxiYn2Pz7ckbj+6V/ioifdn29i2x9+FkbONgj4qPARyHLOgl80IFu89WI2pz8jOOKFl8w5RJ7VURM7Z2Lehvw3a7vAdwhKYB/LVKvn7ObJUp4xXXQ9E8ASLqALNhf2rV5eURslXQicKek+yPih73qKTX9U0T8t5+x23xX1gBdEZJeSPaU67URsT9XdERszX/uBFaTZXbtybnezGpK0inArcAVEfHbru1HSDpq3+/AhcC0I/rdfBlvlmBscoLNTzxWSl0F0j99AjgeuFYSQCe/LVgIrM63tYFbIuL2fudzsJslGDmklTJA11OB9E9vB94+zfaNwJkHHtGbL+PNGsI9u1mKEgff5pqD3SyJoIbz3otwsJulcs9uNv+NTXZ4uKTR+LnmYDdLMHpImyWHH1t1M2bFwW6WzPfsZvNfXZehKcDP2c0awj27Waoh7dkd7GYJxiYmeHjvcC6rVkmwP+OIvbzgxd8/YPtTxxxVQWtsPrnpcycd1PpHWy2WHHnMQT3HweKe3SyRfBlv1gBDPBrvYDdL5bnxZvPf2MQEm/d6uqzZvOcBOrPGGN5XXD2DziyRCn761tM//ZMkfVHShjwF1Fld+1ZIeiDf95Ei7Xawm6UomsK12Ij9DfRO/3QRcFr+WQl8GfanR/9Svv8M4HJJZ/Q7mYPdrCJ5Uoc9PYq8FrgxMncBx0paRLZG/IaI2BgRY8DX87I9+Z7dLMHYRIfNjxUejR80/dNi4OGu75vzbdNtP6dfZQ52swSjrTZLjio8Gj9o+qfpbv2jx/aeHOxmieZwuuxmYGnX9yXAVmB0hu09+Z7dLFV5A3T93AZcmY/Knws8FhHbgJ8Bp0k6VdIocFletif37GYpSpwbXyD90xrgYmAD8CRwVb6vI+ndwPeAFnB9RNzX73wOdrOKFEj/FMC7Zti3huyPQWEOdrME4xMTbHm0oYtXSFoK3AicBEySPV74wqD1mtXRaKuVMhpfK2X07B3gAxGxLs8ZfY+kOyPiNyXUbVY7wzkzvoTR+IjYFhHr8t8fB9aTPfQ3m5/mbjS+VKXes0t6NvAi4O5p9q0km9/LKUuP48zzn1fmqc0A2LBk5jXori3jBDUN5CJKC3ZJRwLfAt4XEQeMYOTTBFcBLDvrlCH977KmG2vyAB2ApBGyQL85Im4to06zOhpttVhy9NFVN2NWBr5nlyTgq8D6iPjc4E0ys4OhjOmyy4ErgJdJujf/XFxCvWb11NQBuoj4McP7NMIsTQzvuvF+EcasITxd1izB+MQEW/Y0eDTerClGWi2WHNPQ0XgzGw7u2c1SDekAnYPdLIGGeDTewW6WysFuNv+NeTTerBlGWy2WHFvOaLykFcAXyNaRuy4iPjtl/4eAN+Vf28DpwAkRsUfSJuBxYALoFFmy2sFulqqEy/iuFE6vJFsy+meSbute9CUirgauzstfCrw/IrozyFwQEbuKntOP3swSlZTYMTWF0+XA1wZpt4PdLEVaYscFktZ2fVZ21TRTaqcDSDqcLAHkt6a05A5J90ypd0a+jDdLMDYxwZbdhQfoeqV/SknhdCnwkymX8MsjYqukE4E7Jd2fJ4qckYPdLMFoq8WS40oZoJsptdN0LmPKJXxEbM1/7pS0muy2oGew+zLerBqFUjhJOgY4D/hO17Yj8pWckXQEcCHw634ndM9ulqqE0fiZUjhJeke+/yt50dcBd0TE77sOXwiszhaJog3cEhG39zung90sVUkz6KZL4dQV5Pu+3wDcMGXbRuDM1PM52M1SeG68WTOMT0ywtfhofK042M0SjLZaLC5nNH7OOdjNUg3pZbwfvZk1hHt2s0QeoDNrgpomgCjCwW6WSDGc0e5gN0s1nLHuYDdL5Xt2s6YY0mD3ozezhnDPbpZgfHyCbTsbPF223yqZZvPFSLvFyccP53TZgS/ju1bJvAg4A7hc0hmD1mtWR0UXmyyw4OScK+OePXWVTLPhVnzByVopI9gLrZIpaeW+VTYf2fVECac1q8a+fG/9PnVTRrAXWiUzIlZFxLKIWHbCgiNLOK1ZBdKWkq6VMgboUlbJNBtq450O23c8VkpdBdI/nU+20OSD+aZbI+IzRY6dThnBvn+VTGAL2SqZbyyhXrPaGWm3WLRg8NH4Iumfcj+KiEtmeezTDHwZHxEdYN8qmeuBb0TEfYPWazbPDTKwPatjS3nOPt0qmWbzVcLg2wJJa7u+r4qIVfnv0w1snzNNHS+R9AuyW+MP5h1p0WOfxjPozFKkDb4Nmv5pHfCsiHhC0sXAt4HTCh57AAe7WYLxzkRZA3R9B7YjYm/X72skXStpQZFjp+NgN0sw0m6x6IRSpsv2HdiWdBKwIyJC0tlkY2y7gUf7HTsdB7tZAlHOhJmC6Z9eD7xTUgd4CrgsIgKY9th+53Swm6Wao/RPEXENcE3RY/txsJulCMBr0Jk1xHDGuoPdLMX4+AQ7tpczXXauOdjNEoyMtFh04jFVN2NWHOxmqXzPbtYANX19tQgHu1miOi5MUYSD3SyVL+PN5r/x8Q7btz1adTNmxcFulmBkpMWihR6NN5v/PEBn1gzZizDDGe0OdrNUwxnrDnazFOPjE2zf6gE6s3lvpN3ipIXDmevNwW6WyJNqzBohhnZSTRnpn8xsCDjYzVJFFPv0IWmFpAckbZD0kWn2v0nSL/PPTyWd2bVvk6RfSbp3ytr0M/JlvFmC8bEJdmz+3cD1FEzh9CBwXkT8TtJFwCqengzigojYVfScDnazBCMjLU5aVMp02f0pnAAk7UvhtD/YI+KnXeXvIlsfftZ8GW+WqpyUzdOlcFrco/zbgO9OacUdku6RtLJIs92zm6UqPhrfK9db4RROki4gC/aXdm1eHhFbJZ0I3Cnp/oj4Ya/GONjNUpWT661QCidJLwSuAy6KiN37mxCxNf+5U9JqstuCnsHuy3izFBGo4KeP/emfJI2SpXC6rbuApFOAW4ErIuK3XduPkHTUvt+BC4Ff9zuhe3azBJ3xCXZs3jNwPQXTP30COB64VhJAJ79SWAiszre1gVsi4vZ+5xwo2CVdDVwKjAH/C1wVEcP5loBZAe2RFgsXHVtKXQXSP70dePs0x20Ezpy6vZ9BL+PvBJ4fES8Efgt8dMD6zOqvpEk1c22gYI+IOyKik38d+DmgWe3ty/XWtGCf4q08/Tng00haKWmtpLWP7HqixNOazaWCgV7DYO97zy7pv4CTptn18Yj4Tl7m40AHuHmmevLni6sAlp11Sv3+J8wK6IxNsOPh3f0L1lDfYI+IV/TaL+ktwCXAy/NE8WbzVnu0xcKTj6u6GbMy6Gj8CuDvySbrP1lOk8xqrMH52a8BDiWbrgdwV0S8Y+BWmdVaA4M9Iv68rIaYDYd6Dr4V4Rl0ZqkmHexm815nrMOOhwqvF1ErDnazBO2RNicuLjgav+7gtiWVg90sie/ZzZrDwW7WAA1+zm7WPA52s/lvfKzD9k2PVN2MWakk2HdvH+Gmf57u3RqzehsZbbFw6TOLFV5/cNuSyj27WYrAk2rMGmNI79m9uqxZoogo9OmnQK43Sfpivv+Xks4qeux03LObJRj/4zg7HtwxcD0Fc71dBJyWf84BvgycU/DYAzjYzRIcffxRXHjleYXK3vzpa3vt7pvrLf9+Y74ozF2SjpW0CHh2gWMPUEmwb9q28YkrP/WGB6o4dwELgDq/6eD2zd5zB61g07aN33vLp/96QcHih/VI/zRdrrfuDK0zlVlc8NgDVNWzP9AjLU6lJK2ta9vA7RtE0TzmvUTEijLaQrFcbzOVKZwnrpsv482qUSTX20xlRgscewCPxptVo2+ut/z7lfmo/LnAYxGxreCxB6iqZ1/Vv0hl6tw2cPsGUZu2Fcz1tga4GNgAPAlc1evYfueUV382awZfxps1hIPdrCEqD3ZJH5QUkoo+uzzoJF0t6f58iuJqSeXk6B2sTcnTI+eKpKWSfiBpvaT7JL236jZNJakl6eeS/rPqtlSl0mCXtJRsyt9DVbZjGrVKRd01PfIi4AzgcklnVNmmKTrAByLidOBc4F01ax/Ae6ndS6dzq+qe/V+AD1OzFBs1TEW9f2plRIwB+6ZH1kJEbIuIdfnvj5MF1eJqW/UnkpYArwauq7otVaos2CW9BtgSEb+oqg0F9UxFPUdmmjZZO5KeDbwIuLvaljzN58k6lcmqG1Klg/qcvVe6Z+BjwIUH8/y9lJWKeo7ManrkXJN0JPAt4H0Rsbfq9gBIugTYGRH3SDq/6vZU6aAG+0zpniW9ADgV+EWeEHIJsE7S2RGx/WC2qV/b9qlZKuoiUysrJWmELNBvjohbq25Pl+XAayRdDBwGHC3ppoh4c8XtmnO1mFQjaROwLCJq8bZUnor6c2SpqCtfXVBSm2yg8OXAFrLpkm8sMmtqLij7i/3vwJ6IeF/V7ZlJ3rN/MCIuqbotVah6gK6urgGOIktFfa+kr1TZmHywcN/0yPXAN+oS6LnlwBXAy/L/r3vzntRqpBY9u5kdfO7ZzRrCwW7WEA52s4ZwsJs1hIPdrCEc7GYN4WA3a4j/B3OReNmnW+ulAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots()\n",
    "\n",
    "with rasterio.open('top-complex.tif') as src:\n",
    "    print(src.profile)\n",
    "    data = src.read(1)\n",
    "    extent = (src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top) #(left, right, bottom, top)\n",
    "    plt.imshow(data.real, vmin=0, vmax=2, alpha=0.5, extent=extent)\n",
    "    \n",
    "with rasterio.open('bottom-complex.tif') as src:\n",
    "    print(src.profile)\n",
    "    data = src.read(1)\n",
    "    extent = (src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top) \n",
    "    plt.imshow(data.real, vmin=0, vmax=2, alpha=0.5, extent=extent)\n",
    "\n",
    "plt.colorbar(ax.get_images()[0])\n",
    "plt.xlim(-4,4)\n",
    "plt.ylim(-3,8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting merged-average-complex.vrt\n"
     ]
    }
   ],
   "source": [
    "%%writefile merged-average-complex.vrt\n",
    "\n",
    "<VRTDataset rasterXSize=\"240\" rasterYSize=\"330\">\n",
    "  <SRS>GEOGCS[\"WGS 84\",DATUM[\"unknown\",SPHEROID[\"WGS84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]]</SRS>\n",
    "  <GeoTransform> -4.0166666666666666e+00,  3.3333333333333333e-02,  0.0000000000000000e+00,  7.9833333333333334e+00,  0.0000000000000000e+00, -3.3333333333333333e-02</GeoTransform>\n",
    "  <VRTRasterBand dataType=\"CFloat32\" band=\"1\" subClass=\"VRTDerivedRasterBand\">\n",
    "  \n",
    "    <PixelFunctionType>average</PixelFunctionType>\n",
    "    <PixelFunctionLanguage>Python</PixelFunctionLanguage>\n",
    "    <PixelFunctionCode><![CDATA[\n",
    "\n",
    "def average(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, buf_radius, gt, **kwargs):\n",
    "    \"\"\" Average overlapping VRTs where (data.reals != 0) \"\"\"\n",
    "    import numpy as np\n",
    "    stack = np.dstack(in_ar)\n",
    "    indValid = (stack.real != 0)\n",
    "    zero2nodata = np.where(indValid, stack, 0)\n",
    "    total = np.sum(zero2nodata, axis=-1)\n",
    "    count = np.sum(indValid.astype(int), axis=-1)\n",
    "    out_ar[:] = total / count\n",
    "    \n",
    "     ]]>\n",
    "    </PixelFunctionCode>\n",
    "    \n",
    "    <SimpleSource>\n",
    "      <SourceFilename relativeToVRT=\"1\">top-complex.tif</SourceFilename>\n",
    "      <SourceBand>1</SourceBand>\n",
    "      <SourceProperties RasterXSize=\"240\" RasterYSize=\"180\" DataType=\"CFloat32\" BlockXSize=\"128\" BlockYSize=\"128\" />\n",
    "      <SrcRect xOff=\"0\" yOff=\"0\" xSize=\"240\" ySize=\"180\" />\n",
    "      <DstRect xOff=\"0\" yOff=\"0\" xSize=\"240\" ySize=\"180\" />\n",
    "    </SimpleSource>\n",
    "    \n",
    "    <SimpleSource>\n",
    "      <SourceFilename relativeToVRT=\"1\">bottom-complex.tif</SourceFilename>\n",
    "      <SourceBand>1</SourceBand>\n",
    "      <SourceProperties RasterXSize=\"240\" RasterYSize=\"180\" DataType=\"CFloat32\" BlockXSize=\"128\" BlockYSize=\"128\" />\n",
    "      <SrcRect xOff=\"0\" yOff=\"0\" xSize=\"240\" ySize=\"180\" />\n",
    "      <DstRect xOff=\"0\" yOff=\"150\" xSize=\"240\" ySize=\"180\" />\n",
    "    </SimpleSource>\n",
    "    \n",
    "  </VRTRasterBand>\n",
    "</VRTDataset>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'driver': 'VRT', 'dtype': 'complex64', 'nodata': None, 'width': 240, 'height': 330, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.03333333333333333, 0.0, -4.016666666666667,\n",
      "       0.0, -0.03333333333333333, 7.983333333333333), 'blockxsize': 128, 'blockysize': 128, 'tiled': True}\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "gdal_vrt_module_0x55cf1cb6f9c0:15: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    }
   ],
   "source": [
    "with rasterio.Env(GDAL_VRT_ENABLE_PYTHON=True):\n",
    "    with rasterio.open('merged-average-complex.vrt') as src:\n",
    "        print(src.profile)\n",
    "        #print(src.read(1)[0:5])\n",
    "        data = src.read(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAEICAYAAACK3Vc9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfbRddX3n8fcnMeADqYAJNISkQY3W4BSkMTpLa7E48tBxgmuNGlqVYXBFW3S0Y1cFO1a0ZaqtD1NXC0xUmvgEphVLtFiljEqdiphQHhIiEoVCTCQ8qCBWJPd+54+973AMN3d/7z37nPPbuZ8Xa6977t77nN/vHvI53332/p3fUURgZmZmZZgz6g6YmZnZo1yYzczMCuLCbGZmVhAXZjMzs4K4MJuZmRXEhdnMzKwgLswdIikkPX3U/TDrCknbJJ04gnYvlvSOYbdrBwb5c8zdISmA5RGxYwhtrQd2RsT/GHRbZtYdw3wdmq38jnlIJD1u1H3IkjR31H0ws8ca5etIl17Dus6FeYAk3SHpbZJuAh6S9DhJR0n6jKR7JN0u6b/17L9K0tcl/VDSbkl/KemgZFtfkfSnkq6T9CNJV0g6vGf730j6fr3tGknH9mxbL+kiSVdKegg4G/ht4A8k/VjS59p7VsyGp87gS+rb59c5+ISkByXdLOkZks6TtEfSXZJe2nPfsyRtr/f9rqTX7/PYf1DndJek1/Veaqoz9Sf17RMl7ZT01rqd3ZLO6nmc35T0L5IeqPtwfs+2ZfXjni3pTuD/SPp7SW/apy83STp9kr9/4v5r637ulvTWnu1TvubU9z1H0m3AbZKuqTfdWL82vGom/19sai7Mg3cG8JvAocA48DngRmAxcBLwFkkn1/uOAb8HLAD+fb39d6fR1muB/wocBewFPtSz7QvAcuAI4Hrgk/vc97eAC4D5wMfq7X8WEYdExMum0Qezkr0M+DhwGPAvwBepXgcXA+8G/nfPvnuA/wj8AnAW8EFJJwBIOgX478BLgKcDv97Q7i8CT67bORv4K0mH1dseosruoVSvFb8zSZH9deBZwMnABuDVExskHVc/7pVTtP9iqvy/FDh34mCF3GvO6cDzgBUR8aJ63XH1a8OnG/5umwEX5sH7UETcFRH/BjwXWBgR746In0XEd4EPA2sAImJLRFwbEXsj4g6qF4mmwPf6eERsjYiHgHcAr5w4LR0Rl0TEgxHxMHA+cJykJ/fc94qI+L8RMR4RP+33jzYr1D9FxBcjYi/wN8BC4D0R8QhwGbBM0qEAEfH3EfGdqHwV+BLwa/XjvBL464jYFhE/Ad7V0O4jwLsj4pGIuBL4MfDMup2vRMTNdfZuAi7lsbk/PyIeql9HrgCWS1peb3sN8OmI+NkU7b+rvv/NwF9TvWHIvub8aUTcX7dtQ+DCPHh39dz+JeCo+rTRDyX9EHg7cCRAfVrt8/Up5weA/0l1JDuTtv4VmAcskDRX0nskfad+3DvqfRbs575mB6q7e27/G3BvRIz1/A5wCICkUyVdK+n+Oqun8WhmjuLnM9OUn/vqg4EJP+lp53mSvlxf3voR8AYem/v///j1wfVG4NWS5lAV2Y83tL/va8NRdduZ1xy/NgyZC/Pg9Q57vwu4PSIO7VnmR8Rp9faLgG9RjXj8BaqirWm0taTn9lKqo/R7qU5Tr6Y67fZkYFm9T+9j7zs838P1bdaSdDDwGeB9wJERcSjVqeKJzOwGju65yxJm7lPAJmBJRDwZuJjH5n7fPG6gGgdyEvCTiPh6Qxv7vjbsqm9nXnP8WjBkLszDdR3wgKoBYU+o38k+W9Jz6+3zgQeAH0v6ZeB3pvn4r5a0QtITqa6X/W39bmA+8DBwH/BEqqPiJncDT51m+2YHioOAg4F7gL2STqW6PjthI3CWpGfVefujPtqaD9wfET+VtIrqQHpKdSEeB95P87tlgHdIemI96PMsYOLa8Exec/zaMGAuzENUF8mXAccDt1O9m/0I1btYgN+nCuWDVNeepzuw4uPAeuD7wOOBiRHfH6M6ffU94Bbg2sRjfRRYUZ9y/7tp9sOs0yLiQar8bAR+QJXLTT3bv0A1uPLLwA5g4h3rwzNo7neBd0t6kKrAb0ze72PAvwM+kdj3q1T9vBp4X0R8qV4/k9ec84EN9WvDK5N9tWnwBCMHCElfAT4RER8ZdV/MZhtJzwK2Agfvcy15kG2+FlgbES+cYp9lVG8C5g2rX9Y/v2M2M5sBSS+XdFD9saf3Ap8bYlF+ItU77XXDaM+Gy4XZbIgkLalH4G5XNY/zmyfZR5I+JGlHPXHECaPoqzV6PdU16O9QfR54umNCZqSe9+Aeqmu9nxpGmza5QeV5YKey6w/g/wUwF/hIRLxnIA2ZdYikRcCiiLhe0nxgC3B6RNzSs89pwJuoPp7zPOAvIuJ5I+kwzrLZ/gwqzwN5x1xPavFXwKnACuAMSSsG0ZZZl0TE7oi4vr79ILCdatamXquBj9UTW1wLHFq/AAyds2y2f4PK86AmJV8F7KhntkLSZXXnbpls5wULFsSyZcsG1BWzyW3ZsuXeiFiY2ffkFz8p7rt/rHG/LTc9vA3onTltXURMeh2wHpjzHOAb+2xazM9P6rCzXrc709eWTSvL4DzbaGTznM0yjC7PgyrMk3Xk5966S1oLrAVYunQpmzdvHlBXzCYn6V+z+953/xjXfXFp435zF93204hYmWj7EKoJLN4SEQ/su3mSu4zq4xONWQbn2UYvm+dslmF0eR7U4K/GjkTEuohYGRErFy5MvWkxG5kAxhP/ZUiaRxXiT0bE5ZPsspOfn6npaB6dqWnYUi8qzrN1RTbLo8zzoApzSS8sZn0LgkdirHFpIklUk7dsj4gP7Ge3TcBr69Gczwd+FBGjOI0NzrIdYLJZHmWeB3Uq+5tU335yDNVsU2tITDNnVrLsEXSDF1B9G9DNkm6o172dav5iIuJiqjmZT6OaqeknVFMojoqzbAeclrIMA8rzQApzROyV9Eaq7zqdC1wSEdsG0ZbZMATBWAsfLYyIr9HwxSRRfYbxnL4ba4GzbAeatrIMg8vzoN4xU3/n6FRf3G3WKeOz9Et2nGU70JSe5YEV5rb8hzmvGHUX7AA1n8N+NbtvAGOFh7kLnGcblGyeu5Dl4guzWSlKP8o2s5zSs+zCbJYQwCP+JjazzutCll2YzRKCKP70l5k160KWXZjNMgLGys6ymWV0IMsuzGYJ1WxBZtZ1XciyC7NZihib+uOKZtYJ5WfZhdksoRowUnaYzaxZF7LswmyWUH32sewwm1mzLmTZhdksabzwo2wzyyk9yy7MZgldOMo2s2ZdyLILs1lCIMYG9i2pZjYsXciyC7NZUumnv8wsp/QsuzCbJQTiZzF31N0wsz51IcsuzGYJ1aQEZZ/+MrNmXciyC7NZUukDRswsp/QsuzCbJUSIsSj7KNvMmnUhyy7MZknjhR9lm1lO6Vl2YTZLqAaMOC5mXdeFLJfdO7NCdGHAiJk160KWXZjNksYK/+yjmeWUnmUXZrOELswWZGbNupBlF2azpPHCR3KaWU7pWe6rMEu6A3gQGAP2RsRKSYcDnwaWAXcAr4yIH/TXTbPRqia+LzvM/XKebTboQpbb6N2LI+L4iFhZ/34ucHVELAeurn8367RAPBJzG5cDgPNsB7RslkeZ50EcNqwGNtS3NwCnD6ANs6GKgLGY07gcgJxnO6BkszzKPPfbcgBfkrRF0tp63ZERsRug/nlEn22YFUCMJ5aOc55tFshleZR57nfw1wsiYpekI4CrJH0re8c6+GsBli5d2mc3zAYr4EB9R9zLebYDXhey3FfvImJX/XMP8FlgFXC3pEUA9c89+7nvuohYGRErFy5c2E83zIZijDmNS5c5zzZbZLI8yjzPuGVJT5I0f+I28FJgK7AJOLPe7Uzgin47aTZqgRiP5qWrnGebLbJZHmWe+zmVfSTwWUkTj/OpiPgHSd8ENko6G7gTeEX/3TQbrQAeKXx+3T45zzYrdCHLM+5dRHwXOG6S9fcBJ/XTKbPyqPjvcO2H82yzR/lZLvuwwawQQfmzBZlZsy5k2YXZLKn0o2wzyyk9y2UfNpgVIkKMx5zGJUPSJZL2SNq6n+1PlvQ5STdK2ibprFb/GLNZLJvlTJ4HlWUXZrOEasBIa1P4rQdOmWL7OcAtEXEccCLwfkkH9dN/M6tks5zM83oGkGWfyjZLUWuTEkTENZKWTbULMF/VEOlDgPuBva00bjbrlZ9lF2azhGrASOq61AJJm3t+XxcR66bZ3F9SfX54FzAfeFVEjE/zMcxsEtPIMvSf5xll2YXZLCk5E9C9Pd/MNFMnAzcAvwE8jWp6zH+KiAf6fFwzY1pf+9hvnmeUZV9jNksY8sxfZwGXR2UHcDvwy209uNlsNuSZv2aUZRdms6Rx5jQuLbmTelIPSUcCzwS+29aDm812mSy3lOcZZdmnss0SIuCR8XYKr6RLqUZoLpC0E3gnMK9qJy4G/hhYL+lmQMDbIuLeVho3m+W6kGUXZrOE6vRXayM5z2jYvovqSyTMrGVdyLILs1lS6bMFmVlO6Vl2YTZLmOZHLMysUF3IsguzWUp7p7/MbJTKz7ILs1nSeOGnv8wsp/QsuzCbJVQjOdNzYZtZobqQZRdms4SJSQnMrNu6kGUXZrOk0k9/mVlO6Vl2YTZL6MJITjNr1oUsuzCbJZU+ktPMckrPsguzWUKE2Ft4mM2sWRey7MJsllT66S8zyyk9yy7MZglduC5lZs26kGUXZrOk0sNsZjmlZ7nxRLukSyTtkbS1Z93hkq6SdFv987CebedJ2iHpVkknD6rjZsOU/XL10jnPNttlszzKPGeugK8HTtln3bnA1RGxHLi6/h1JK4A1wLH1fS6UVPYUK2ZJ46hx6YD1OM82y2WyPMo8NxbmiLgGuH+f1auBDfXtDcDpPesvi4iHI+J2YAewqqW+mo1MBOwdn9O4lM55ttkum+VR5nmm15iPjIjdABGxW9IR9frFwLU9++2s15l1XhdOVc+Q82yzSulZbnvw12R/bUy6o7QWWAuwdOnSlrth1q4uzK87AM6zHXC6kOWZvle/W9IigPrnnnr9TmBJz35HA7sme4CIWBcRKyNi5cKFC2fYDbPhiVDj0lHOs80qmSyPMs8zLcybgDPr22cCV/SsXyPpYEnHAMuB6/rrolkZSh4s0ifn2WaV0gd/NZ7KlnQpcCKwQNJO4J3Ae4CNks4G7gReARAR2yRtBG4B9gLnRMTYgPpuNjQR5V+XynCebbbrQpYbC3NEnLGfTSftZ/8LgAv66ZRZecRYB0ZdN3GezcrPsmf+Mkvq8DVkM+tRepZdmM0SujC/rpk160KWXZjNMqK6NmVmHdeBLLswmyV1eNS1mfUoPcsuzGYJ0YEBI2bWrAtZdmE2Syr99JeZ5ZSeZRdms6TSR3KaWU7pWXZhNkuIKD/MZtasC1l2YTZLKv0jFmaWU3qWXZjNkkq/LmVmOaVn2YXZLCEQ44WP5DSzZl3IsguzWVLhB9lmllR6lss+bDArRbT3fcySLpG0R9LWKfY5UdINkrZJ+mprf4fZbJfMcibPg8qyC7NZViSWnPXAKfvbKOlQ4ELgP0XEsdRfw2hmLclkOZfn9Qwgyz6VbZbU1kcsIuIaScum2OW3gMsj4s56/z2tNGxmQPlZ9jtms4QAxsfVuAALJG3uWdbOoLlnAIdJ+oqkLZJe2+ofYzaLZbPcUp5nlGW/YzbLCCB3lH1vRKzss7XHAb8KnAQ8Afi6pGsj4tt9Pq6Z5bMM/ed5Rll2YTZLGuJnH3dSvSA8BDwk6RrgOMCF2awFpWfZp7LNstob/NXkCuDXJD1O0hOB5wHbW3t0s9muvcFfTWaUZb9jNkvJfxyq8ZGkS4ETqa5f7QTeCcwDiIiLI2K7pH8AbgLGgY9ExH4/jmFm01F+ll2YzbJaekccEWck9vlz4M/badHMfk7hWXZhNssIiPGyJ743s4QOZNmF2Syt7DCbWVbZWXZhNssqfYJdM8spPMuNo7InmwtU0vmSvlfP/3mDpNN6tp0naYekWyWdPKiOmw3d8EZlD4zzbMYwR2XPSObjUuuZfC7QD0bE8fVyJYCkFcAa4Nj6PhdKmttWZ81GZmJSgqalfOtxnm02y2Z5hHluLMwRcQ1wf/LxVgOXRcTDEXE7sANY1Uf/zIoR0byUznk2y2V5lHnu5xrzG+t5PzcDb42IHwCLgWt79tlZr5vaI1sZ//7ySTft+OAb+uiiGTz9965t3imj8JGcfXKerRNayXPhWZ7pzF8XAU8Djgd2A++v10/210563CFp7cTE4PfcNzbDbpgNj6J56Sjn2WaVTJZHmecZFeaIuDsixiJiHPgwj57e2gks6dn1aGDXfh5jXUSsjIiVC5/iy1ZWuMIHi/TDebZZJZvlrhVmSYt6fn05MDHCcxOwRtLBko4BlgPX9ddFsxKUPVikH86zzS7JLI8wz43XmPczF+iJko6nOqa4A3g9QERsk7QRuAXYC5wTET6vZQeGjr4j7uU8m1F8lhsL837mAv3oFPtfAFzQT6fMijQ+6g70z3k2o/gse+Yvs4zpfbm6mZWqA1l2YTZL6vCoazPrUXqWXZjNsgoPs5klFZ7lmX6O2czMzAbA75jNkko//WVmOaVn2YXZLCMofho/M0voQJZdmM2yCj/KNrOkwrPswmyWVPrpLzPLKT3LLsxmWYWH2cySCs+yC7NZVuFhNrOkwrPswmyWMOqvgTOzdnQhyy7MZlmFj+Q0s6TCs+zCbJZU+lG2meWUnmUXZrOswsNsZkmFZ7n4wvydV1086i5Yx815822Trpe0Jf0gHbgu1QXOs/Wr7zx3IMvFF2azYhQeZjNLKjzLLsxmSSr8y9XNLKf0LPvbpczMzArid8xmWYWf/jKzpMKz7MJsltGBASNmltCBLLswm2UVHmYzSyo8yy7MZlmFh9nMkgrPsguzWYIofySnmTXrQpY9KtssIx6d/H6qJUPSJZL2SNrasN9zJY1J+s9t/AlmRjrLmTwPKsuNhVnSEklflrRd0jZJb67XHy7pKkm31T8P67nPeZJ2SLpV0smZjpgVLxJLznrglKl2kDQXeC/wxRn1dfLHdJbNIJflXJ7XM4AsZ94x7wXeGhHPAp4PnCNpBXAucHVELAeurn+n3rYGOLbu8IV1x8y6raXCHBHXAPc37PYm4DPAnpl1dlLOshm0VpgHleXGwhwRuyPi+vr2g8B2YDGwGthQ77YBOL2+vRq4LCIejojbgR3AqmyHzEqVPPW1QNLmnmXttNuRFgMvB1qdWNpZNqtM41R2X3meaZanNfhL0jLgOcA3gCMjYjdUgZd0RL3bYuDanrvtrNft+1hrgbUASxd7DJp1QO4d8b0RsbLPlv4X8LaIGJMG872xbWa5fjzn2bojf9mp3zzPKMvpBEk6hOrt+Fsi4oEpGplsw2OehohYB6wDWHnc4wsfvG6zXgx1JOdK4LI6YwuA0yTtjYi/a+PB284yOM/WIR3IcqowS5pHFeRPRsTl9eq7JS2qj7AX8ej5853Akp67Hw3syv8dZoUaUrmJiGMmbktaD3y+xaLsLJsVnuXMqGwBHwW2R8QHejZtAs6sb58JXNGzfo2kgyUdAywHrsv8EWYla/HjUpcCXweeKWmnpLMlvUHSGwbaf2fZDGj141IDyXLmHfMLgNcAN0u6oV73duA9wEZJZwN3Aq8AiIhtkjYCt1CNAj0nIsb66aRZEVo6yo6IM6ax739pp1XAWTarFJ7lxsIcEV9j8mtNACft5z4XABdkO2FWvOl9TrlIzrIZnciyh0+aJYjyv5HGzJp1IcsuzGZJpYfZzHJKz7ILs1lW4WE2s6TCs+zCbJZVeJjNLKnwLLswm2VM4+NQZlawDmTZhdksq/Awm1lS4Vl2YTZLKv3L1c0sp/QsuzCbJZV++svMckrPsguzWUYHJiUws4QOZNmF2Syr8DCbWVLhWXZhNkvowmxBZtasC1l2YTZL0njhaTazlNKz7MJsltGB61JmltCBLLswmyWVfvrLzHJKz7ILs1lW4WE2s6TCs+zCbJZU+lG2meWUnmUXZrOswsNsZkmFZ9mF2Swjyp/Gz8wSOpBlF2azhC589tHMmnUhyy7MZllReJrNLKfwLLswmyWVfpRtZjmlZ9mF2SyjA5MSmFlCB7LswmyWVPqAETPLKT3LLsxmSaWH2cxySs/ynKYdJC2R9GVJ2yVtk/Tmev35kr4n6YZ6Oa3nPudJ2iHpVkknD/IPMBuKoBow0rQUzFk2I5/lEeY58455L/DWiLhe0nxgi6Sr6m0fjIj39e4saQWwBjgWOAr4R0nPiIixNjtuNmylDxhJcJbNKD/Lje+YI2J3RFxf334Q2A4snuIuq4HLIuLhiLgd2AGsaqOzZiMViaVgzrJZLZPlEea5sTD3krQMeA7wjXrVGyXdJOkSSYfV6xYDd/XcbSeThF/SWkmbJW2+5z4fgFvZJiYlaFq6os0s14/nPFsnZLM8yjynC7OkQ4DPAG+JiAeAi4CnAccDu4H3T+w6yd0f8ydGxLqIWBkRKxc+Ze60O242VBFovHnpgrazDM6zdUgyy6PMc6owS5pHFeRPRsTlABFxd0SMRcQ48GEePcW1E1jSc/ejgV3tddlsRAo+9ZXlLJvR/VPZkgR8FNgeER/oWb+oZ7eXA1vr25uANZIOlnQMsBy4rr0um41Gyae+Mpxls0rpp7Izo7JfALwGuFnSDfW6twNnSDqe6rjiDuD1ABGxTdJG4BaqUaDneBSndV4AHTlVPQVn2awDWW4szBHxNSa/1nTlFPe5ALigj36ZlafsLDdyls1qhWd5WqOyzWaztk591SOf90jaup/tv12PkL5J0j9LOq7Nv8NstmvrVPagsuzCbJbU4ijO9cApU2y/Hfj1iPgV4I+Bdf313Mx6tTgqez0DyLLnyjbLaHGUZkRcU3+OeH/b/7nn12upRkObWRs6kGUXZrOEalKCVJoXSNrc8/u6iOjnHe/ZwBf6uL+Z9ZhGlqHdPKez7MJslpX7Rpp7I2JlG81JejFVmF/YxuOZWS3/7VKt5Hm6WXZhNkuaxlF2/21JvwJ8BDg1Iu4bWsNms0DpWXZhNssY4kxAkpYClwOviYhvD6dVs1miA1l2YTZLaW/uXEmXAidSXb/aCbwTmAcQERcDfwQ8BbiwmqyLvW2dHjez8rPswmyW1dLpr4g4o2H764DXtdKYmT1W4Vl2YTbLCFB+wIiZlaoDWXZhNssa4oARMxugwrPswmyWVXaWzSyr8Cy7MJslabzw819mllJ6ll2YzTKC6UxKYGal6kCWXZjNEkQMdVICMxuMLmTZhdksq/Awm1lS4Vl2YTbLKjzMZpZUeJZdmM0yOnBdyswSOpBlF2azpNJHcppZTulZdmE2S4niT3+ZWUb5WXZhNssIig+zmSV0IMsuzGZZZZ/9MrOswrPswmyWVPpnH80sp/QsuzCbZRUeZjNLKjzLc5p2kPR4SddJulHSNknvqtcfLukqSbfVPw/ruc95knZIulXSyYP8A8yGIgLGxpuXgjnLZuSzPMI8NxZm4GHgNyLiOOB44BRJzwfOBa6OiOXA1fXvSFoBrAGOBU4BLpQ0dxCdNxuqiOalbM6yGeSyPMI8NxbmqPy4/nVevQSwGthQr98AnF7fXg1cFhEPR8TtwA5gVau9NhuFgoOc4Syb1bpemAEkzZV0A7AHuCoivgEcGRG7AeqfR9S7Lwbu6rn7znrdvo+5VtJmSZvvuW+sn7/BbPACGI/mpXCDyHL9uM6zdUM2yyPMc6owR8RYRBwPHA2skvTsKXbXZA8xyWOui4iVEbFy4VN8dsxKFxDjzUvhBpHl+nGdZ+uIZJZHmOdpjcqOiB9K+grV9aa7JS2KiN2SFlEdgUN1VL2k525HA7va6KzZyATFD+6aDmfZZq0OZDkzKnuhpEPr208AXgJ8C9gEnFnvdiZwRX17E7BG0sGSjgGWA9e13XGzoSv4mlSGs2xWK/wac+Yd8yJgQz0acw6wMSI+L+nrwEZJZwN3Aq8AiIhtkjYCtwB7gXMiwhedrPsKL7wJzrIZFJ/lxsIcETcBz5lk/X3ASfu5zwXABX33zqwY5b8jbuIsm0EXsuyZv8wyAij8q+LMLKEDWXZhNssq/CjbzJIKz7ILs1lKFD+S08wyys+yC7NZRkB04HPKZtagA1l2YTbL6sDMXmaWUHiWXZjNsgq/LmVmSYVn2YXZLCOi+JGcZpbQgSyXUZjnPZs5v7h51L0wm1rhR9nFcJ6tdIVnuYzCbFa8IMY86ZVZ95WfZRdms4yJr4ozs27rQJZdmM2yCv+IhZklFZ7l1Pcxm812AcR4NC4Zki6RtEfS1v1sl6QPSdoh6SZJJ7T5t5jNZtksZ/I8qCy7MJtlRPLL1XPWU30P8v6cSvUVi8uBtcBFffXdzB6VzXIuz+sZQJZ9Ktssqa0BIxFxjaRlU+yyGvhYRARwraRDJS2KiN2tdMBslis9y0UU5i1bttwr6SHg3hF3ZYH7MPL2h9mHX8ru+CA/+OI/xt8uSOz6eEm9nxVaFxHrptmvxcBdPb/vrNd1ojBv2bLlx5JuHXE3ZtO/Y/ehksrzNLIM/ed5RlkuojBHxEJJmyNi5Sj74T6Mvv1S+rCviJjqdFXbNFkXhth+v24d9f+/Ev4NuQ/l9KFXF7Lsa8xm5dkJLOn5/Whg14j6YmYzN6MsuzCblWcT8Np6ROfzgR/5+rJZJ80oy0Wcyq5N9zrcILgPo28fyujDwEi6FDgRWCBpJ/BOYB5ARFwMXAmcBuwAfgKcNZqezlgJ///ch4r7MECDyrKi8DlDzczMZhOfyjYzMyuIC7OZmVlBRl6YJZ0i6dZ6yrJzh9juHZJulnTDxOfUJB0u6SpJt9U/D2u5zcdM3zZVm5LOq5+XWyWdPMA+nC/pe/VzcYOk0wbchyWSvixpu6Rtkt5crx/qc2Htc55nV56d5QGJiJEtwFzgO8BTgYOAG4EVQ2r7DmDBPuv+DDi3vn0u8N6W23wRcAKwtalNYEX9fBwMHFM/T3MH1Ifzgd+fZN9B9WERcEJ9ez7w7bqtoV/VzisAAAIrSURBVD4XXtpdnOfZl2dneTDLqN8xrwJ2RMR3I+JnwGVUU5iNympgQ317A3B6mw8eEdcA9yfbXA1cFhEPR8TtVKP6Vg2oD/szqD7sjojr69sPAtupZsMZ6nNhrXOeZ1meneXBGHVh3t90ZcMQwJckbZG0tl53ZNSfMat/HjGEfuyvzWE/N29U9e0nl/Scdhp4H1TNM/sc4BuU81zYzDjP5fwbHnqeneX2jLowj3LqwRdExAlU3/5xjqQXDandrGE+NxcBTwOOp5rD9f3D6IOkQ4DPAG+JiAem2nWQ/bDWOM/7d0Dn2Vlu16gL88imHoyIXfXPPcBnqU6n3C1pEUD9c88QurK/Nof23ETE3RExFhHjwId59NTSwPogaR5VkD8ZEZfXq0f+XFhfnOcC/g0PO8/OcvtGXZi/CSyXdIykg4A1VFOYDZSkJ0maP3EbeCmwtW77zHq3M4ErBt2XKdrcBKyRdLCkY6i+z/O6QXRgIkC1l1M9FwPrgyQBHwW2R8QHejaN/LmwvjjPBfwbHmaeneUBGfXoM6rpyr5NNTrvD4fU5lOpRgbeCGybaBd4CnA1cFv98/CW272U6tTSI1RHjmdP1Sbwh/Xzcitw6gD78HHgZuAmquAsGnAfXkh1+uom4IZ6OW3Yz4WX9hfneXbl2VkezOIpOc3MzAoy6lPZZmZm1sOF2czMrCAuzGZmZgVxYTYzMyuIC7OZmVlBXJjNzMwK4sJsZmZWkP8HseCZQhB9HJEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,(ax1,ax2) = plt.subplots(1,2, figsize=(8,4))\n",
    "plt.sca(ax1)\n",
    "plt.imshow(data.real, vmin=1, vmax=2)\n",
    "plt.colorbar()\n",
    "plt.title('real part')\n",
    "\n",
    "plt.sca(ax2)\n",
    "plt.imshow(data.imag, vmin=1, vmax=2)\n",
    "plt.colorbar()\n",
    "plt.title('imaginary part');"
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}