notebooks/isce2_usage.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a8392f58-d0ff-4f7d-8e4c-dc17a65e5bcf",
   "metadata": {},
   "source": [
    "# Using sardem for ISCE2 processing\n",
    "\n",
    "While [ISCE2](https://github.com/isce-framework/isce2) has the `dem.py` script to create DEMs, you can use the DEMs produced by `sardem` in ISCE by setting the `-isce` or `--make-isce-xml` (for example, if you'd like to use the Copernicus DEM's, which aren't currently available through the ISCE DEM script). \n",
    "\n",
    "The `-isce` flag will produce the necessary  \".xml\"/\".vrt\" files used by ISCE.\n",
    "\n",
    "Note that ISCE is *not* installed as a dependency of `sardem`; you must separately `conda install -c conda-forge isce2` or see the ISCE README for full installation instructions to use this flag."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15e1d0a6-7fe6-4937-97e4-357daf661333",
   "metadata": {},
   "source": [
    "## Differences between sardem and ISCE's `dem.py` script:\n",
    "The main differences between how a user requests a DEM with `sardem` and with `dem.py` are due to `sardem` following GDAL's bounding box convention.\n",
    "\n",
    "\n",
    "1. From the help string, one usage of `dem.py` is :\n",
    "\n",
    "> Stitch the requested files and apply EGM96 -> WGS84 correction (-c)\n",
    ">\n",
    "> dem.py -a stitch -b 31 33 -114 -113 -r -s 1 -c\n",
    "\n",
    "We see that the order of the boundaries for ISCE is \"SNWE\" (south, north, west, east), while `sardem` uses `--bbox left bottom right top`, which means the order  is (west, south, east, north).\n",
    "The `dem.py` script also currently allows only integer bounds, so that the entire SRTM tiles are stitched together. `sardem` allows any floating-point bbox.\n",
    "\n",
    "\n",
    "2. In GDAL's bounding box convention, the four numbers refer to the **outer edge** of the image. In ISCE, the XML files point to the **center** of the pixels.\n",
    "\n",
    "What this means in practice is that to completely replicate the output of \n",
    "\n",
    "`dem.py -a stitch -b 31 33 -114 -113 -r -s 1 -c`\n",
    "\n",
    "\n",
    "the bounding box needs to be shifted up/left 1/2 a pixel:\n",
    "\n",
    "`sardem --bbox -114.000138888888 31.000138888888 -113.000138888888 33.000138888888`\n",
    "\n",
    "This is because the SRTM lat/lon numbers in the tile names refer to the pixel centers, so the tiles actually extend out 1/2 a pixel in all directions.\n",
    "Also, the `dem.py` script leaves off the bottom/right row of pixels in the tile (hence a shift up+left, rather than a padding in all directions).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "04e83f19-7268-4a8e-af2e-90cb99f090f2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This is the Open Source version of ISCE.\n",
      "Some of the workflows depend on a separate licensed package.\n",
      "To obtain the licensed package, please make a request for ISCE\n",
      "through the website: https://download.jpl.nasa.gov/ops/request/index.cfm.\n",
      "Alternatively, if you are a member, or can become a member of WinSAR\n",
      "you may be able to obtain access to a version of the licensed sofware at\n",
      "https://winsar.unavco.org/software/isce\n",
      "curl -n  -L -c $HOME/.earthdatacookie -b $HOME/.earthdatacookie -k -f -O https://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11/N32W114.SRTMGL1.hgt.zip\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100 8288k  100 8288k    0     0   486k      0  0:00:17  0:00:17 --:--:-- 2401k--:--     0\n",
      "curl -n  -L -c $HOME/.earthdatacookie -b $HOME/.earthdatacookie -k -f -O https://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11/N31W114.SRTMGL1.hgt.zip\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100 5846k  100 5846k    0     0   297k      0  0:00:19  0:00:19 --:--:-- 1352k\n",
      "API open (R): ./demLat_N31_N33_Lon_W114_W113.dem\n",
      "API close:  ./demLat_N31_N33_Lon_W114_W113.dem\n",
      "Writing geotrans to VRT for ./demLat_N31_N33_Lon_W114_W113.dem\n",
      "GDAL open (R): ./demLat_N31_N33_Lon_W114_W113.dem.vrt\n",
      "API open (WR): demLat_N31_N33_Lon_W114_W113.dem.wgs84\n",
      "\n",
      "    << Geoid Correction I2 SRTM>>\n",
      "\n",
      "Jet Propulsion Laboratory - Radar Science and Engineering\n",
      "\n",
      "\n",
      "Sampling Geoid at grid points -  Longitude Samples:    13 Latitude Lines:    23\n",
      "Corner Geoid Heights (m) =  -34.35 -35.82 -32.83 -31.68\n",
      "\n",
      "Correcting data to geoid height...\n",
      "\n",
      "At line:      512\n",
      "At line:     1024\n",
      "At line:     1536\n",
      "At line:     2048\n",
      "At line:     2560\n",
      "At line:     3072\n",
      "At line:     3584\n",
      "At line:     4096\n",
      "At line:     4608\n",
      "At line:     5120\n",
      "At line:     5632\n",
      "At line:     6144\n",
      "At line:     6656\n",
      "At line:     7168\n",
      "GDAL close: ./demLat_N31_N33_Lon_W114_W113.dem.vrt\n",
      "API close:  demLat_N31_N33_Lon_W114_W113.dem.wgs84\n",
      "API open (R): demLat_N31_N33_Lon_W114_W113.dem.wgs84\n",
      "API close:  demLat_N31_N33_Lon_W114_W113.dem.wgs84\n",
      "Writing geotrans to VRT for demLat_N31_N33_Lon_W114_W113.dem.wgs84\n",
      "Writing geotrans to VRT for /Users/staniewi/repos/sardem/notebooks/demLat_N31_N33_Lon_W114_W113.dem.wgs84\n",
      "N32W114.SRTMGL1.hgt.zip = succeded\n",
      "N31W114.SRTMGL1.hgt.zip = succeded\n"
     ]
    }
   ],
   "source": [
    "!dem.py -a stitch -b 31 33 -114 -113 -r -s 1 -c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "bde05925-30aa-4847-84f3-a9abf3a24ff5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[08/17 12:34:01] [INFO dem.py] Bounds: -114.000138888888 31.000138888888 -113.000138888888 33.000138888888\n",
      "[08/17 12:34:01] [INFO download.py] /Users/staniewi/.cache/sardem/N32W114.hgt already exists, skipping.\n",
      "[08/17 12:34:01] [INFO download.py] /Users/staniewi/.cache/sardem/N31W114.hgt already exists, skipping.\n",
      "[08/17 12:34:01] [INFO dem.py] Cropping stitched DEM to boundaries\n",
      "[08/17 12:34:03] [INFO dem.py] Rate = 1: No upsampling to do\n",
      "[08/17 12:34:03] [INFO dem.py] Writing DEM to sardem_version.dem\n",
      "[08/17 12:34:03] [INFO dem.py] Writing .dem.rsc file to sardem_version.dem.rsc\n",
      "[08/17 12:34:03] [INFO dem.py] Creating ISCE2 XML file\n",
      "This is the Open Source version of ISCE.\n",
      "Some of the workflows depend on a separate licensed package.\n",
      "To obtain the licensed package, please make a request for ISCE\n",
      "through the website: https://download.jpl.nasa.gov/ops/request/index.cfm.\n",
      "Alternatively, if you are a member, or can become a member of WinSAR\n",
      "you may be able to obtain access to a version of the licensed sofware at\n",
      "https://winsar.unavco.org/software/isce\n",
      "Writing geotrans to VRT for /Users/staniewi/repos/sardem/notebooks/sardem_version.dem\n",
      "Creating output file that is 3600P x 7200L.\n",
      "Processing egm_sardem_version.dem [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.\n"
     ]
    }
   ],
   "source": [
    "!sardem -isce --bbox -114.000138888888 31.000138888888 -113.000138888888 33.000138888888 -o sardem_version.dem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "9cd35689-3b10-4bfe-9e86-66f56cb0b762",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sardem_version.dem     sardem_version.dem.vrt\n",
      "sardem_version.dem.rsc sardem_version.dem.xml\n"
     ]
    }
   ],
   "source": [
    "!ls sardem_version*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1036f8e9-def6-494d-abf7-e85aa5179b57",
   "metadata": {},
   "source": [
    "The \".xml\" and \".vrt\" files were created since the `-isce` flag was passed.\n",
    "\n",
    "We can see from the \"startingvalue\"s that `sardem` has shifted the convention to be pixel-centeres for these ISCE files, and that they match the `dem.py` values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "81782a95-9f47-49b1-b962-329239d7fede",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30-        <property name=\"size\">\n",
      "31-            <value>3600</value>\n",
      "32-            <doc>Coordinate size.</doc>\n",
      "33-        </property>\n",
      "34:        <property name=\"startingvalue\">\n",
      "35-            <value>-114.0</value>\n",
      "36-            <doc>Starting value of the coordinate.</doc>\n",
      "37-        </property>\n",
      "38-    </component>\n",
      "--\n",
      "59-        <property name=\"size\">\n",
      "60-            <value>7200</value>\n",
      "61-            <doc>Coordinate size.</doc>\n",
      "62-        </property>\n",
      "63:        <property name=\"startingvalue\">\n",
      "64-            <value>33.0</value>\n",
      "65-            <doc>Starting value of the coordinate.</doc>\n",
      "66-        </property>\n",
      "67-    </component>\n"
     ]
    }
   ],
   "source": [
    "!grep -n4 'starting' sardem_version.dem.xml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "6c0a5b92-1fdc-48f6-8077-7e14e1ea35b9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "33-        <property name=\"size\">\n",
      "34-            <value>3600</value>\n",
      "35-            <doc>Coordinate size.</doc>\n",
      "36-        </property>\n",
      "37:        <property name=\"startingvalue\">\n",
      "38-            <value>-114.0</value>\n",
      "39-            <doc>Starting value of the coordinate.</doc>\n",
      "40-        </property>\n",
      "41-    </component>\n",
      "--\n",
      "62-        <property name=\"size\">\n",
      "63-            <value>7200</value>\n",
      "64-            <doc>Coordinate size.</doc>\n",
      "65-        </property>\n",
      "66:        <property name=\"startingvalue\">\n",
      "67-            <value>33.0</value>\n",
      "68-            <doc>Starting value of the coordinate.</doc>\n",
      "69-        </property>\n",
      "70-    </component>\n"
     ]
    }
   ],
   "source": [
    "!grep -n4 'starting' demLat_N31_N33_Lon_W114_W113.dem.wgs84.xml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "b068433a-3ab2-4ec9-bbb3-88fc70e2f7be",
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "\n",
    "ds = gdal.Open(\"demLat_N31_N33_Lon_W114_W113.dem.wgs84.vrt\")\n",
    "dem_isce = ds.GetRasterBand(1).ReadAsArray()\n",
    "ds = None\n",
    "\n",
    "ds = gdal.Open(\"sardem_version.dem\")\n",
    "dem_sardem = ds.GetRasterBand(1).ReadAsArray()\n",
    "ds = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "2b0bb74a-ad51-4580-b99e-9b9f8ea2ec40",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 0.012789197530864197)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(dem_sardem - dem_isce).max(), abs(dem_sardem - dem_isce).mean(),"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "54e8bf0b-aee6-4b22-a116-9f58a7aa2bc6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWEAAAHwCAYAAAB3xd4xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABMbklEQVR4nO3deXCc+XnY+e/z9oWDIO6DOGZ4gQeI4THkjEZrr9fRYY2dxKPdKiXjrKNJoqqpaJVE3kpsz0RVdlK7s6Wss06sdaStWVvr0Ua2osRWacq2ZI/H9jpOZM3wJghe4IAkQJA4CRBH4+juZ/94X3CaJG5299tvv8+n6i10/7pf9O8Fuh/88Duen6gqxhhj/OH4XQFjjAkzC8LGGOMjC8LGGOMjC8LGGOMjC8LGGOMjC8LGGOMjC8LmISKiIrK3COrxP4rIH/tdD2PyzYKwKUqq+k1V/Yl8fX8R+dcick1EpkXksoh8Nl+vZcxaon5XwPhDRKKqmvK7Hj6aBf4mcBV4Dvi+iPSp6n/1t1ombKwlHEAi8osicttrxV0RkY975c+LyA9EZFJE7ojIr4tIPOs8FZEviMg14JpX9vPec4dE5B888joJr8V4S0SGReT/EpFy77EfF5FBEfkFERnxvsenReSnROSqiEyIyD9/gmv8eyLyl95tEZF/473OlIicF5Fu77FyEfk/ROSm99hfZtXxBRH5r97P45yI/Pjy91fVX1bVy6qaUdUfAv8Z+OhW62vMVlkQDhgR2Q/8I+A5Va0CPgXc8B5OA/8z0IAbUD4O/E+PfItPAx8BukTkReCfAZ8EOoFPPPLcfwXsA44Ce4E24JeyHm8ByrLK/2/gZ4HjwH8L/JKI7H6Cy132E8CPeXWpAf42MO499q+91/tvgDrgF4CMiLQBfwD8r175PwN+V0QaH/3mXtB+DriYg7oaszmqakeADtxgOIIbMGPrPPfngO9k3VfgY1n3vw58Oev+Pu85ewHB/Zd9T9bjHwX6vds/DiSBiHe/yjv3I1nPPwV8eovX+feAv/Rufwy32+AFwMl6juPV4cgK5/8i8P8+UvZHwCsrPPct4PuA+P37tSN8h7WEA0ZV+3CD678ARkTkWyLSCiAi+0Tk90XkrojcB/433FZxtoGs262P3L+ZdbsRqABOef/OT+IGquyW5Liqpr3bSe/rcNbjSWDbo9cgIk+JyMzysYFr/lPg14F/BwyLyJsist27tjLg+gqnPQ18ZrnuXv1/FNjxSF1+BegG/paqWjYrU3AWhANIVX9bVX8UN9AobrcBwNeAy0Cnqm4H/jlui/ah07Nu3wE6su4/lXV7DDeIHlLVGu+oVtXHguoW6n9LVbctHxs85yuqehw4hNti/3mvjvPAnhVOGcBtCddkHZWq+uXlJ4jIvwR+EvgJVb3/pNdlzFZYEA4YEdkvIh8TkQRuAEri9gWD2yVwH5gRkQPA59f5dt8G/p6IdIlIBfDLyw+oaga3j/ffiEiT99ptIvKp3F7R+kTkORH5iIjEcLtI5oG0V8evA78qIq0iEhGRj3o/m38P/E0R+ZRXXuYNJrZ73/N14O8An1TV8VVe2pi8syAcPAngy7itwLtAE26LF9zBp78DTOMG0P+w1jdS1e8B/xb4U6DP+5rtF73yv/K6N/4E2J+Li9ik7bjXcw+3y2Qcd0AO3Gu+ALwPTOD+V+Co6gDwEu7PZhS3ZfzzfPie/99wW/7XsrpGtjybw5itEusGM8YY/1hL2BhjfFTwICwiL3oLDPpE5LVCv74xxgB44wTveQt5LnoDtYWvRyG7I0Qkgjvf85PAIG4/3s+oam/BKmGMMbgrMYFKVZ3xBn3/Eviiqv5VIetR6Jbw80Cfqn6gqovAt3AHT4wxpqDUtTxPPeYdBR8kK3QCnzYeXhwwiLuE9iEi8irwKkBlZdnxAweeevQpwaMZWFqASBQiMb9rY8yG3Lhxl7GxqUfnmufU9pouTS3N5vz7JuduXcSdzrjsTVV9M/s53n/np3BXif47dfOIFFShg/BKv8zH/vJ4P6g3AU6c2K/vv//mYycFUwamRyGdhuomEEtiZ4rbc8+9mvfXSC3NcqD7F3P+fc+894V5VT2x1nO8FZ9HRaQG+I6IdKtqT84rs4ZCR4FBHl6h1Q4MFbgOPnKgqhHSKVhMQqISm6Biwk4ExMlrY3tdqjopIn8OvAgUNAgXOgK8D3SKyC4vxeLLwNsFroPPHIjEHwTgxUwllsLDmMITkUavBbycSe8TuMv+C6qgLWFVTYnIP8LNZhUBvq6qIU0f6AbemJNkKVOOqkM8MouQ8blexhSeRHxpCe8A3vL6hR3g26r6+4WuRME7JVX1D4E/LPTrFishQ9yZRXFYTFeSJkbCuU9EwrzphTH5p6rngWN+18NGhoqEkCERmUZxmE9XkdEYich9orLod9WMyS8BnPB2yVkQLjJChvLIFIpDMl1DMhOnLDJFzEmuf7IxgSR+dUcUBQvCRUrIUBGZQCMOc6k6ZlJNVEQmSESm/a6aMSaHLAgXOSFDZXSMCiaYWWri/tIOKmNjVEQm/K6aMblRBFPU/BTejpiAETJUxe7SUNbHYrqSu8luZlJNflfLGPOErCUcMEKGmvgAab3D1GIHAwvPUZu4ybboiN9VM2bLwtwnbC3hgIpIirpEP60VZ0imarg9d4xkutrvahljNslawgEXkRSNZVdZypQzsbCLe7qThrI+4k7uE6IYkxc2Rc2UgpiTpLm8l8VMJWPze4nIEnWJfpvaZoqe2BQ1U0riziytFedIpqsZTnYhkqEh0WdT24wpUhaES1R5ZIr2ylPMpJq4NXOQsqjSUHaV8siU31Uz5mEhn6JmQbjEbYuOsLd6jMnFDvom97ItHqW54prNMzamSFgQDgEhQ238JtvrbjO+sJdL4zupKdtLS8UHVEbH/K6eMaHuEw7vkGQIRSRFU9llDjdcIhFJc260nWtTzzObavC7asaElrWEQyjmJGmvPEV9WTWDM/s5NdzGjm272VF5wxZ9mMKzPmETVuWRKTqr36OpopUbUzs4ebeFHduepn1bv3VTmMIRgUh4/ym3IGyojg3RXT/C+MJe+u45DM+6LePWSuszNibfLAgb4MP+4trmcu4mu+mfXOD2dCttVbtp29ZnsylMXll3hDGemJOko/J96hIN3Li/k4H789yZ6aCtai+tlVdsnrExOWZB2KyoMjpGV90EzZW7uH6viptT89ye3snOmjKay3ptBZ7JGZFwT1GzIGxWJWRoSFynuqmc4WQXH0ym6JtIciuym53VUZrKL1uiIGOekAVhs67lKW21XhfF8OwiVyeWuBHZy64ah6ayy5YoyDwR6xM2ZgOWuyiaKvbywWQFs0sZroynuRndz+6aFA1lV213aLN5Eu4sauGdnGe2RMjQWHaVZ5uusKumHEdgPpWhd8zhzMghRuf3kVb7227MRtmnxWxJzEmyq+qHNJS10DfZyr35FNOLaS6MllFbdpid1ePUxAcQMn5X1QRBiJO6h/fKTU5Uxe5yuOE8B+qFuPcv5b35FGeGq+mdeNY2IzVmHdYSNk8sIilaK85Rk6jj5v093JlZAGB4dpHRuRY6tu+0OcZmTWEemLOWsMmZisgEB2pP8UzjPGVR962VUbg5Nc+p4Z3cnjtGSuM+19KY4mItYZNTywN31c2V3Jw+xO3peTIKi2nlyniakbJu9tXespwU5kMhX6xhLWGTF3Fnls7q9zjSNEVl7MO32b35FKeG2xiZP+Bj7YwpHhaETV7Vxm9yrOkaHdvLHpSlMkrPaJzr95+37gnj7rbs5P4ICgvCJu+WW8WHm2YfzKAAt6/43Gg3c+k6H2tnfCe4+YRzfQREcGpqAq8hcZ3jzTdorvyw9Tu1kOL0cAcTC7t8rJkx/rEgbAqqPDJFV91pDtQLUe9fxsW0cn60isHZ46i9JUPJuiOMKSAhQ2vFOY41DVMVjwDuVLarE0v0TZ2wZc8mVCwIG99Uxe5ytPHyQ4N2A/fn6Z04wlKm3MeamYLypqjl+ggKC8LGVzEnyd7qk3Q3Lj7onhidW+L82EEW0lU+186Y/LMgbHwnZGgqu8yzzXcedE9MLaQ4O7qH2VSDz7Uz+WdT1IwpCtuiIxxtvExbVQKA2aUM50bbmV5q8blmJq+sO8KY4hFzkuyrOcWBenmQq/jcaDOTix1+V82YvLBhaFN0lmdPVEQ76B2vZz6V4fxoHc80OtTGb/pdPZNjAojlEzam+NTEBzjW1E9tWZRURjk3Um2LOkzJsSBsilp5ZIpnGnro2F5GRuH8aJUF4lIj4DiS8yMoLAiboheVRfZWn+RAvfvBskBsSokFYRMIy/3EhxuncUQsEJeU3LeCrSVsTJ7UJfp5tvkO5VGH86NVNmuiBIhNUTMmWNz5xNepTkQ5P1rH1FKr31UyZsssCJtASkSmeaahh8aKOBdGG21X54Cz7ghjAigqi+yrOUNbVRnnR1stObwJJFusYQItIil2Vr1PPHKEi2NPc7hhiURk2u9qmc3wpqiFlbWETeAJGdoqztBaBb0TnZaP2ASKBWFTMlrKL5BRGJo75ndVzCYI4e4TtiaDKRkRSXGg7hanhtuoTTSxLTrid5XMhgiO5Y4wpjRURsfYVRPnykS7dUuYQFg3CIvI10VkRER6ssrqROQdEbnmfa3Neux1EekTkSsi8qms8uMicsF77CsiEpz/F0ygtFacIZXJcDf5jN9VMRsgAk5Ecn4ExUZawr8FvPhI2WvAu6raCbzr3UdEuoCXgUPeOV8VkYh3zteAV4FO73j0exqTExFJsb9unOv3MiTT1X5Xx5g1rRuEVfUvgIlHil8C3vJuvwV8Oqv8W6q6oKr9QB/wvIjsALar6g9UVYFvZJ1jTM7VxAdorIjTN7kXtV63ohfmgbmtvjubVfUOgPd1eblSGzCQ9bxBr6zNu/1o+YpE5FUROSkiJ0dHp7ZYRRN2e6p7uTefYtwS/Zgilusmwkp/fnSN8hWp6puqekJVTzQ22r+TZmvizix7a+Hy+DZSGve7OmY1gm30uQXDXhcD3tfluUCDQHZaq3ZgyCtvX6HcmLxqLr9IIuLY3OEiJpbKckveBl7xbr8CfDer/GURSYjILtwBuPe8LotpEXnBmxXx2axzjMmbiKTorB3jxuQ8C+kqv6tjzGPWnUgpIr8D/DjQICKDwC8DXwa+LSKfA24BnwFQ1Ysi8m2gF0gBX1DVtPetPo8706Ic+J53GJN3NfEB6subuTVzkM7q9/yujnlUyHNHrBuEVfVnVnno46s8/w3gjRXKTwLdm6qdMTmyq/oDTt7toG1bHRWRRyf7GOMfW1JkQqEiMkHH9k76p3ZzqM6CcDERCNTiilyzCZQmNNorz3NvfskSwJuiYi1hExoxJ8muGof+qQ6eqbfkPkUj5H3C1hI2odJSfoHpxbS1houKm0Ut18e6ryrSISJ/JiKXROSiiHyxABf7GAvCJlQikmJv7TzXJ22XZkMK+KeqehB4AfiCl/+moCwIm9BpLLvKfCrN9FKL31UxuANzfqyYU9U7qnrauz0NXGKNdAr5YkHYhI6QYU/tNDfuF/zzZgqrYTkHjXe8utoTRWQncAz4YcFq57GBORNK9Yl++iePMJOyHTh8l7+BuTFVPbHuy4tsA34X+DlVvZ+PiqzFWsImlIQMu2ru0z9lfcNhJiIx3AD8TVX9PT/qYC1hE1rLreHZVAOV0TG/qxNqfizW8PLY/CZwSVV/teAV8FhL2ITWcmt4aHa331UJNRHfkrr/CPB3gY+JyFnv+Kn8Xu3jrCVsQq0+0c/NqWqWMuXEnKTf1TEFpKp/ycq5zgvKgrAJNSHDU9tnGF/YQ0t5z/onmDwIVv7fXLPuCBN6dYkPGJ5N2F50xhfWEjahF5EUzZUL3F9qoTpmG74UmljuCGNMQ9lVRuba13+iMTlmLWFjgKgsUhFbYjFTSdyZ9bs6oSMbSLhTqiwIG+NpLLvKxMIuG6ArNLGk7sYYIO7MksrYAJ0pLGsJG5Olrqyf2VSD5ZMoILEpasaYZRWRCWaWLOG7KRxrCRvziPLoBCmNE5VFv6sSDjZFzRiTbXvsLjNLzX5Xw4SEtYSNeYSQ8bsKoRPmlrAFYWNWsC02bHOGC0QAR8IbhK07wpgVRGWRpUy539UwIWAtYWNWYaktC8QG5owxK4k7s9YaNnlnLWFj1iCS9rsKJU8QItYSNsasxOYKm3yzlrAx68pg7ZX8CnOfsAVhY9ajGW8nMgvE+SBiU9SMMWsRxw3ExuSBtYSN2TDrlsgX644wxqzBse4IkzcWhI3ZEAfFsbwSeRLmlrD9aTdmg9JqbRaTe/auMmaDIpIirVEikvK7KiXF3fLe71r4x4KwMRskZEhlKolEpv2uSokRIjZFrYipgvXDGWNKVPG3hEW8yfLF//fClL54ZJaFdBUJaw3njGADc8VPHNuG3BQFIcN8psrvapgSUvwtYQBvalBK4wgZGxgxvoo58zZAl0uWTzg4orJIRmMsZir9rooJsfLIJMl0nd/VMCUiIC3hD8WcJGmN2v5fxjdChpjYrhu5YnvMBVBEUsSdWVIa97sqJqQSkWlIzWMzd8yTClxLOFtUFrOWklpyFVNIGUhOQVWj3xUpCWHuEw50EAY+XMtvOV9NQTmQqITJu1DTgr3vnoBIqINw6bxzxIF0CjSF/YtoCiJegU4Mue87Y7Yo8C3hDzkQiUJq0Q3IkSil9DfGFCMH2d4Ao/3Q0om937ZGwDb6LB0ORONu18TCLNYiNnlX04LeugFL837XxARUiQVheBCIowmYvWfdEya/onGoqUZvX8beZ1vniOT8CIoSDMLgdk3EoaIaZsbdLgr7gJi8cJDmp+HWICQtn4TZvHWDsIh0iMificglEbkoIl/0yutE5B0RueZ9rc0653UR6RORKyLyqazy4yJywXvsKyJ5/nMlUXcKUXLK+4BYIDZ5sL0RolG0/wL2Hts88ZYt5/oIio20hFPAP1XVg8ALwBdEpAt4DXhXVTuBd737eI+9DBwCXgS+KiIR73t9DXgV6PSOF3N4Latw3ECcScHUMPYhMTknDtLeCiNj9h7bIgvCa1DVO6p62rs9DVwC2oCXgLe8p70FfNq7/RLwLVVdUNV+oA94XkR2ANtV9QeqqsA3ss7JMwcqa6FsG0wMQtq6J0wuOdC8BxwHvd7rDgwbs0Gb6hMWkZ3AMeCHQLOq3gE3UANN3tPagIGs0wa9sjbv9qPlK73OqyJyUkROjo5ObaaKa/Am11e3wL0hWJzDArHJmUQlPNUO96dh7Bb23to4kdy3gkuqJbxMRLYBvwv8nKreX+upK5TpGuWPF6q+qaonVPVEY2P1Rqu4Ad6AXcNTMHcf5qawD4vJDQdp3ee1hq96g8HGrG9DQVhEYrgB+Juq+nte8bDXxYD3dcQrHwQ6sk5vB4a88vYVyn3guEtNrZ/Y5FK8wm0Nz8/D3WvY+2rjbIraGrwZDL8JXFLVX8166G3gFe/2K8B3s8pfFpGEiOzCHYB7z+uymBaRF7zv+dmsc3zgwLb6D/uJ1ZaemiflIDs63dbwjZtel5cxa9vIsuUfAf4ucEFEznpl/xz4MvBtEfkccAv4DICqXhSRbwO9uDMrvqCqae+8zwO/BZQD3/MOH3n9xJEYjA+6reNonJKdPm3yL1EJrS0wOITe6kX2nsDeT2sL+x5z6wZhVf1LVu7PBfj4Kue8AbyxQvlJoHszFcw/B6JlUN/uZsSqqHY/SPbBMVviIB0H0aG7MHQXdtxzZ+bY+2l1YrkjDLgLO2pbYX7GXWVn/Xlmq5Zbw5kMer3H79qYImdB+CEOVDeDE3VbxRaIzZa4rWEcBybuuVMi7b20quXtjWxgzngct0uiYrsN2JmtS1RCizt1Xq9fsgUcZlUWhFfkuNONtje5A3a2l5jZtKzW8MwsjN7A3kOrscUaZkVZA3b3R2yFndm88ipoagBA+2wBh1mZBeH1SBTq2mF20lbYmU1ykKe81vDiIty5gr1/ViDujyjXR1AEqKp+ctyZE0vztsLObE5FNTTUAaA3B7wdX0w2G5gzG+TNnIhE3QE7C8RmQxzk6YPuzVQKHbiEvXdMNgvCm+ItdU5UwsgHNnPCbExlLdR5ex4M3XW33bJA/BBHcn8EhQXhTfNyE1c1wp1rlpvYbICD7Drg3rQFHOYRFoS3xHFHvhs63MEW28POrKeqEWq8tKwT92wxUBbrEzZb5M0lbt7jBmKbwmbWITv3P7itfRdtAYcBLAg/IQdiZbBjP9ztsylsZg0OVDfBtkr3ri3geIi1hM0T8BZ1tHWhI/2W/MesThxkd1Zr+IM+W8AB7jxhG5gzTywSRzq60bEB6+8zq/Dmm1eUu3dtBw6DBeHc8gIx90dss0ezskdbw7YDBwJERHJ+BIUF4VyLxKGjG+an3bnEFojNQxx3o1kvpwSpFDp4GXufFJ6IfF1ERkTE1zmDFoTzQaLQfggyaft306zAQfYcgai3sc3gkDeoG14+9Qn/FvBiXi9sAywI540DLZ3u/nW3L9nqOvOwRCXS6XVLZDLojYvYH+vCUtW/ACb8rocF4bxyoHGnm8Rl8KIFYpPFgaZdHy7gGBmD6VHCGIiF3E9P86aoNYjIyazjVb+vdSUb2W3ZPBFvRNyJwEAPtHW5/cbGiIPsO4qe/M9ua/haD3LsvwMJX9soT/N6x1T1RD6+cS6F77ftCy8DW10bOtBj+SaMx91KS/bsce/en7b96ELIgnDBuBnYpGmXG4gt34QBHowdeCvptO8SpMPVbSW2WMMUjtfyad6DDvZaIDauSBTZfxgcB+aSMHYTe1/kn4j8DvADYL+IDIrI5/yoh/UJF5ybgU12dKK3epCOLjf/hP09DDHHzbL2VDvcuIVev4bUd7jL4UPCj1wPqvozBX/RFdgn3xeOO0Wp/YA7jzhk/36alThIexeUldl+dCFjQdg3XirMlk73A2fT10w0juzrArz96EK0nNn6hI1PvFSYTbtg+LoF4tBz3J29mxrc5cy3eglDa9iSuhufeS3iujYYtw1ETdaS5qG7oV/OHAYWhIuCF4ir6m2eqHHHC/Z2hmg5c95WzAWCBeGi4Q7WUVFtieFDz3G3zdpeFerlzGFhQbioeIE4VgYLs9gHL8TEQfYfBcdBr/f6XZu8ssUapsh4gRhsMUeoOVBZi+zeBZNTMGHjBaXKgnBR8gJxJmU78oaa424iu73KXc5cwu8F6xM2RcgbrEunsBZQiEWibrfE/EJJ785s3RGmSDkQtbSX4fZht4TtzlyaLAgXPQdwSKkF4/DyuiXiMXdRT4m1hm2xhgmEqCyykK7yuxrGL163hA4sL2curUAcZhaEAyQembVAHFpet8SefejlUyU3cyYiuT+CwoJwgAgZ4pFZljLlflfF+MLds1Dq6tCrJy37XomwIBwwQoaos2CBOLQcaDsI0Sh6/XRJJH1yF2tYn7AJkOVAbIN1ISVRZPdRmEui/WdLIhBbEDaBI2SIyiJptc1RQikaRw6ecHNLDIYhyU/psk9wwEUkheIg9iEMGW93lu7j6NkfIpE4tO4niO0qd4qa37XwT/B+Y+YxQga1X2UIeTt4H+hG+67ByAdYizh4rCVcIqwlHFYONDyF7JlFL/cisTKobSVY7atg9eHmWpB+U2ZdGW+QxgJyuDhuV0RrC3rxrOUfDhgLwqVGln+l9iEMFYkiu47Ctkr0wilvW6SAvAcsn7ApHW6eCTRj2dfCKBJFup5zE8H3vBeYjQEsd4QpPcutYQvEIeOmP5Xu47C4hF58D5bmsfdAcbMgXJIciETdFrF9CEPmwxkTzMyiV04FYHmzImRyfgSFBeGS5eUitkAcQsszJvbA2HjJLG8uVTZFraR5u3PMTbl3Y2XY392w8GZMzE7D0F0o64GnDlOsv3+R8DYS1v2NiEiZiLwnIudE5KKI/EuvvE5E3hGRa97X2qxzXheRPhG5IiKfyio/LiIXvMe+IhKg3vPAcqCiGuanrUUcNss5Jmqq0Q9uwNAV7PdffDbyZ3EB+JiqHgGOAi+KyAvAa8C7qtoJvOvdR0S6gJeBQ8CLwFdFJOJ9r68BrwKd3vFi7i7FrM6BqkaYmSi5PLRmHcs5JsrK3FV1RbhPnaA4pHN+BMW6QVhdM97dmHco8BLwllf+FvBp7/ZLwLdUdUFV+4E+4HkR2QFsV9UfqKoC38g6x+Sd466kujdk+5SFipdj4tAxAPRSD0zepegCsWRyfgTFhjqIRCQiImeBEeAdVf0h0KyqdwC8r03e09uAgazTB72yNu/2o+WmYByoa4fxAUhbIA4P9z8hOdAFmQzaexZm71FsgTisNhSEVTWtqkeBdtxWbfcaT1+pn1fXKH/8G4i8KiInReTk6OjURqpoNioSh/oOd8NIGzEPEQeadiO7d8LiInrxZBEt5rApahumqpPAn+P25Q57XQx4X0e8pw0CHVmntQNDXnn7CuUrvc6bqnpCVU80NlZvpopmI6JxqG93B2osEIeIAx3d0NTgJoS/dNK6porARmZHNIpIjXe7HPgEcBl4G3jFe9orwHe9228DL4tIQkR24Q7Aved1WUyLyAverIjPZp1jCsrtJ6SmxUbMw0aiSOezsL0KJqfcveqK4A+x9QmvbQfwZyJyHngft0/494EvA58UkWvAJ737qOpF4NtAL/B94AuqujxU+XngN3AH664D38vhtZhNcXfvpaIa7l7DAnGIxMqQQ8/DtkoYGUNvnsd+//5Zd7GGqp4Hjq1QPg58fJVz3gDeWKH8JLBWf7IpKG/GxNgtNyF4026KdTK/ySVvxsThj6LnfwA3bkG83LedOZanqIWVfeJCz13iCrjB2FpEIeEl+zn8UTf9Zd81+/37xIKwYXnknMVkUc4hNfnySCC+3ONPQnixPmFjeJBrYGbcPSwQh8Ry+suPuIH44hlITlPo379NUTMGAAfaDsLE7WDtzGCekBeID30EyhLu1LWleb8rFRoWhM3DJAptXehIfxFN5jf553izJj7i7sxxrZDpL9W6I4x5SCSOtHeht69Y5rVQyQrEi0ve1DWTbxaEzcqicaSjC719OQA7M5jcyQrEk1MF65YQ0jk/gsKCsFmF92Hc0Qm3e4tiVZUplKxAnLLfe77ZzhpmDd7y5voOuH0J2g9hf7fDwg3ElFXk/ZUExQlQH26uWRA26/CWN6cW3eXNLZ1YIA4L58Odu/MsSFPKcs0+TWYDHKhudpe2FuHODMYEmQVhs0FeQnjN2Ko6k2M2Rc2YDfKWN8/es1V1xuSI9QmbTfKWNw9edHfpKK/C/pabJyEQqCllKxGRZzfwtCVVvfBooQVhs3kShdaD6EAP0nbAHUW3QGzC7f/Dzbe+0jZuy3YBOx8ttCBstiYSRdoOoLcvIx3dbqvYmC3RQPXhruJ9Vf3YWk8QkT9dqdyCsNkib0J/8x4Y8uYQi72dzNY4AR9fWC8Ar/Uc+x/SPAHH7ROubYU7tkWSMSLyIyJS6d3+WRH5VRF5eq1zLAibJ+TAtnooq7Q5xGbLSiif8NeAORE5AvwCcBP4xlonWBA2OeDNIc6kbQ6xCbuUqirwEvBrqvprQNVaJ1gnnskRB5r3uFPXYgl3qbP9jTcbIFISA3PLpkXkdeBngR8TkQgQW+sE+5SYHFremWMQFuewFrEJob8NLACfU9W7QBvwK2udYC1hk1veHGLuXIEd+23qmtmQoC/WAPBavf9eVT+xXKaqt1inT9iCsMm9SNzNtjZ8HXZ02tQ1s47SSGWpqmkRmRORalWd2uh59ukw+RGNQ0MHjN2Cxp1Yz5cJiXnggoi8A8wuF6rqP1ntBAvCJk/cHXzZ3ujOmKhpwQKxWU0J5RP+A+/YMAvCJo+8nTnSKTfr2rZ6LBCbUqaqb4lIOfCUql7ZyDn2iTB55kBFtbtDw8IsNmPCPEpKKJ+wiPxN4Czwfe/+URF5e61zLAibAvC2SFqad7dJskBsSte/AJ4HJgFU9Sxu9rRVWRA2BeItb16YdXfnMCaLX1vei8iLInJFRPpE5LUcXEpqhZkRutYJ1idsCshrESenLRm8yeLPijlvXu+/Az4JDALvi8jbqtr7BN+2R0T+DhARkU7gnwD/da0T7FNgCszLvGYr6oz/ngf6VPUDVV0EvoWb8+FJ/GPgEO6qud8GpoAvrnWCtYSND9xcxCzOudPYrC0QepKfLqoGETmZdf9NVX0z634bMJB1fxD4yBO+5l9X1S8BX1ouEJHPAP9xtRMsCBt/SNRdWTc35c6esEBscm9MVU+s8fhKWxGt2X+7Aa/zeMBdqewBC8LGP5Gom3Ft9p5lXQs7fwZrB4GOrPvtwNBWvpGI/CTwU0CbiHwl66HtQGqtc4v/Xa8ZrO+wVHndEtE4TI9iv2dTYO8DnSKyS0TiwMvAmnN61zAEnMRdtnwq63gb+NRaJxZ/S1gc0upm4orImn9QTCB5q+pSi24grmokCG0Dk0vqS0tYVVMi8o+APwIiwNdV9eIWv9c54JyI/DZuXC2tFXMRSZHKlDObavC7KiYvvKlr6TRMDWMt4pBRQDX3x0ZeWvUPVXWfqu5R1TdycDUvUqor5hKRaUQyjM7vQ4NTbbNhjpvkZzHp5pmwQGyC6V/w+Iq5nWudEKhoVhGZYFtsmIHZ51hIr7ltkwkkx017OTNhgThsNJP7wx8rrZhbU6CCMEB5ZIqW8h4+uH+AycWO9U8wAePtVTd51501YYHYBMtDK+ZE5P+kFFfMxZ1Z9lb3cHummcHZ46S1+McXzSZIFHbsR0dvukucLRCXOIVMJveHP7JXzP0OcB/4ubVOCGQQBog5SQ7UnmV2yaFn/Ahz6Tq/q2RyKRJH2rvQgUuWAtMEhqrOqeqXVPU5VT3h3Z5f65zABmFwZ03sqzlFbXmU08MdTCysmTHOBE00juzsRvsvWK6JUlcifcIickJEfk9ETovI+eVjrXMC/3+8kKGj8n0icoSese3srT3CjooLpbRdSoi5WyTJ04fQD84je466izuC3XYwj1J/5gnnyTeBnwcusMFWQ0m8m4UMbRVnOFC/QN89uDH9nPUTlww365p07Eevn7Wk8KbYjarq26rar6o3l4+1TiipSNVUdplY49NcHKthbuko+2ouEnOSflfLPDF3MYe07UE/OIvsedZN/mNKR+m0hH9ZRH4DeBd3cA4AVf291U4oqSAMUBu/ydGmJOdHWzk7eoBDDTepiEz4XS3zxByoakRSS2i/F4il5N6+Jvj+PnAAiPHhv2wKhCcIA2yLjnCsaYHe8T2cHu6gu6GSmvjA+ieaIudAbSuSWkRvnkeePmyBuFT4N6Us146o6jObOaEk+oRXUh6Z4nDDJWrLYpwdqWU42eV3lUxOuKvqpKwSBi9i/cOmyPyViGwq2JRsEIYP5xJ3bC/j0niU/umP2IBdSXCgpRPEgaErWCAOujwk79lgAp88+FHgrLd56HkRuVDyU9TWE5EUu7efpCx6hGsTSWYWj7C/9jJxZ9bvqpkn4kDbQbdbYuQDaNpNibcpSpdSSgNzL272hFC8a5ensB1qmOfefIpzo/uYSTX5XS3zpCSKPH0YnRyDsVtYi9j4RUROA2RPS3t0itrycx5V8i3hbI1lV4k3tXJxrInTwzvoaqiiIXHd72qZJyFRZM+z6JX3kGjcTYcZjrZFaQl+S/jgOt0OAlSv9MCGg7CIRHC377itqn9DROqA/4CbK/MG8LdU9Z733NeBzwFp4J+o6h955ceB3wLKgT8Evqha2M6b6tgQx5pm6R3fQ89oJXtrj9NWecZW2AVZJIp0Poteeh+JxmBbPRaITYEd2MBz0isVbuad+kXgUtb914B3VbUTd2LyawDeyODLuJmEXgS+6gVwgK8BrwKd3rHp/pNcWJ450VgR5+rEEpfvHWcpU+5HVUxOuHvVyf5n0b4ey7wWOIpqJudHQa9glW6IR47Blc7dUBAWkXbgrwO/kVX8EvCWd/st4NNZ5d9S1QVV7Qf6gOdFZAewXVV/4LV+v5F1TsEtz5zYW1fOnZkFzo4esO2TAs3dq046D6NXTlvCHxMYG20J/1vgF3j4Xd2sqncAvK/LI11tQPbKiEGvrM27/Wj5Y0TkVRE5KSInR0c3laR+UyKS4qnKH3KoIcXsUppTw22MLezJ2+uZfPOWN+86gF49Y3kmgkIppXzCm7ZuEBaRvwGMqOqpDX5PWaFM1yh/vFD1TS8X54nGxhX7snOqubyXo033iDpCz2glA7OWACi4HKhuRlo60L7TkLYduouflkwqy63YSEv4R4CfFpEbwLeAj4nIvweGvS4GvK8j3vMHgex9h9qBIa+8fYXyolATH+BYUz9V8SjXJhbonTjCYqbS72qZLXGg4Slkey168zyoBWJTvNYNwqr6uqq2q+pO3AG3P1XVnwXeBl7xnvYK8F3v9tvAyyKSEJFduANw73ldFtMi8oKICPDZrHOKwvKAXVtVgtG5Jc6MdDK11Op3tcyWONC6H4lGbXlzEJTOirlNe5J5PF8GPiki14BPevdR1YvAt4Fe4PvAF1R1eWrG53EH9/qA68D3nuD18yLmJNlXc4p9dTGSqQxnhhu4m+z2u1pmSxzo6EaTSbh7DQvEphhtquNTVf8c+HPv9jjw8VWe9wbwxgrlJ4Gij2hChvbKU5RF93B5fBu9Yw6LdR+ho/J9m08cNBJFdh9FL7+PxMuhrh2bQ1yEAtSHm2v2blxDQ+I6x5puUxWP0DeR5OrkcVJqycQDJxpH9h9Hb12H6VGsRWyKiQXhdVRGxzjaeJkd2xLcnl7g0sQztrAjcLzFHAdPoNd7YW4KC8TFxGZHmHW4CzvcfuLx5BLnxw6STOd/6pzJJW8xx76jaN95W8xRTGyesNmI5X7iw43TJFNpzozsskxsgeNARbW7mKPvrLeYwxh/WRDepLpEP882D5CIOJwd2cHkYsf6J5ki4u1V17obvWFziIuDdUeYTaqITHCksYfmygTnR+sYnd/nd5XMpjhQ04LUtdgcYuM7C8JbFJVF9lafZF9dmkvj5QzNHUHtxxkgjjtdrawKhq9jgdhn1hI2WyFkaCnv4UjTKAP3hRvTlnMiWNxNQ9EM3BvCArHxgwXhHHATxV9jbinN1cljNpc4ULxNQ+emYGYcC8Q+sWXL5knFnVm66k6zPZHi4vgzLKSr/K6S2TA3zwSTd2FhFgvEBaZqU9RMbixvKPr09hEu39tnSeKDRKLQdhBGb9jUNVNQFoTzoCY+wIHaqwzO7OLe4tN+V8dslERhx34Y6bepa4VmA3Mm1xKRafZWnyGZqmFkfiN7AJqiEIlC49NuILZuCVMAFoTzKCIpWivOEXdmuZvstilsgeDmmaCmBSYGsUBcINYSNvlUEx+gLtHPSPKAzZwIBDfPBGVVlnXN5J0F4QKJO7M0lV9memmHzZwIBDfPBOJAchoLxHkU8tkRtrKggIQMtfGbJNPVLKSrSESm/a6SWZMD2+rd1nAqAdE41m7JkwAFzVyzd5QPyiNTRJ2k5SUOBDfhD8mpQPUzmuCwlrBPIpIiIinSGiUiNh2quHmBeGbcbRlb2yXH1FrCxj8RSVmLOBAcqKy1XTlMzllLuAjEnCTJdDVlkWnbSLSYSdSdNbE4B/EKrA2TIwpkgpPrIdfsXVQkyiNTzKXqLAtbsYtEIRL3ljbbH0zz5OwTX0Qqo2NML7VQER2zfuKi5biBOJ1yB+rE2jE5YX3CplhUxe4yvbTDFnUUNS8QZzJYa9g8KWsJF6Ga+ABjC3vYHrtL3Jn1uzpmRV4g1gyId99sjYZ7doQF4SJVn+hndH4f1bHbtqijaDluALZA/ORCHITtXVOkhAyNZVcZX9jNXLrO7+qYVTluv/DSPNY1YbbCWsJFTMiwo+ICg7PH0YRDZXTM7yqZFXmZ15LTUF6FtW22wFrCplgJGVorznB3bjfTSy1+V8esyoGySpi9h7WIzWZYEA6AiKTYWXWaodkO26mjmEnUTX85NYwF4k0IeRY1C8IBEZEUe6rPMTTTyNjCHr+rY1YTibopMO8NYYHYbIQF4QCJyiL7a89zZ6aG4WSX39UxK/L6hytrLBBvRkZzfwSEBeGAicoiB2p7GZop4/bcMb+rY1bkuLklEpUwdgsLxOtQrDvCBEvMSdJdf4mRWeHW7Eds77qi5O3MkSi3verMmuzTG1AxJ8mh+iuMzi7RN3XCEv8UJS8PMVggXpMNzJmAijuzPNNwhcn5JXonjlhe4qLkQF07LC1YH7FZkQXhgIs7sxxpvMp8KkPP+EHbRLQoOdC8x51DPHkXC8SPU83k/AgKC8IlIO7M0t3wAYvpDGdH99gy56LkQNtBdGLI3TjUArHxWBAuEeWRKQ439gNweriDqaVWn2tkHiNR5Klu9PZ1W1mXzRZrmFKxHIhjjnBupJHJxQ6/q2QeFY0juw+j13vcXBMWiF0WhE2pWA7EiYhwfrTOAnHRcecQy55u9Mppd786C8ShZkG4BFkgLnbuzs2ycz96+ZSlwYSiWzEnIp8RkYsikhGREzm6yhVZEC5RFoiLnQM1LUhzK3rllLtnnSkmPcD/APxFvl/IgnAJs0Bc7Nypa7K9Cr1+GjSkgbgIB+ZU9ZKqXsnRFa7JgnCJs0Bc7Bzo6IZUCr15ntB3S+RWg4iczDpe9btCK7G1riFQHpniaON1zo/t5fxoHYcb3c1ETZGQKLLvBHr+v0DsErQdJHTto/zMZhhT1VX7c0XkT4CVdkr4kqp+Nx8VWokF4ZBIRKY50niVc6P7LBAXo2gcOXgCPfcDJJaApt2EJhAvZ1Er9MuqfqLgL7qCkPyWDXy4xLk86ljXRNFxoLwKOXAYvXrZ8kyEiAXhkFkOxFXxiAXiouPNmNh3AL10HmbGCUcgLr6BORH570VkEPgo8Aci8kc5udQVWBAOobgzS3f9JQvERcmBpt1IRxvacwoWZglHIC4uqvodVW1X1YSqNqvqp/L1WhaEQ2o5MXx1ImqBuOg40H4I6mrQi++FYzFHkbWEC8mCcIi5ieEvUFtmgbjoSBTZ8yxEo7aYo8RZEA65qCzSVXeO+vKYBeJiE4kiB5+D+fnSnkNchIs1CmlDQVhEbojIBRE5KyInvbI6EXlHRK55X2uznv+6iPSJyBUR+VRW+XHv+/SJyFdERHJ/SWazIpLiQO1ZmivjnB+tY2Jhl99VMsDyzs1y6HkYGYOhK5RsIC6y3BGFtJmW8F9T1aNZk59fA95V1U7gXe8+ItIFvAwcAl4EvioiEe+crwGvAp3e8eKTX4LJhYik2Ft9hraqBD1j2y0QFw1v6lr3cfRGv+3eXIKepDviJeAt7/ZbwKezyr+lqguq2g/0Ac+LyA5gu6r+QFUV+EbWOaYIRCTF7u0n6dheZoG4qDiwrR45eBi92lt6U9dsy/sNUeCPReRU1vrrZlW9A+B9bfLK24DspViDXlmbd/vRclNEhAw7q95nV02cnrHtjC3s8btKBgAHaluRPZ02da3EbHTZ8o+o6pCINAHviMjlNZ67Uj+vrlH++DdwA/2rAE891bzBKppcETJ0VL6PI8foHdvGwfp9NJZd9btaZjnr2uIceukk0v1RiJb5Xakc0EC1XHNtQy1hVR3yvo4A3wGeB4a9Lga8ryPe0weB7CH2dmDIK29foXyl13tTVU+o6onGxuqNX43JqbaKM+yrS3NpvJyR+QN+V8cAD+YQV5SjH5wNb/rLErJuEBaRShGpWr4N/ARuwuO3gVe8p70CLGcdeht4WUQSIrILdwDuPa/LYlpEXvBmRXw26xxTpFrKezhYn+TyeILhZJff1THgziHefRTmknD7EiXRLRHiPuGNdEc0A9/xZpNFgd9W1e+LyPvAt0Xkc8At4DMAqnpRRL4N9AIp4Auqmva+1+eB3wLKge95hylyjWVXiTTsond8O2k9QmvFOb+rZKJxpOs59Ox/QcqqoOEpAjvtX0HTwZlSlmvrBmFV/QA4skL5OPDxVc55A3hjhfKTQPfmq2n8Vpfo55nGVi6MNgIWiP3nbRh66AR64X2kvAoqawlsIA4x+42ZDauODXG06Q4D98W6JoqCt2Hovi700ulg55iwxRrGbMy26AjPNN7i1v0Eo/P7/K6OwYG6dmSHt2GoBjQIh5gFYbNpFZEJDjf0MTRTZQs6ioIDrfuhLAEDPQSuNawK6TwcAWFB2GxJIjJNV10vw3MNlvSnGEgU2XUUnZyC0RsELhCHmAVhs2UxJ0lnzTnG53cwvbTSfommoCJRZP+z6M3rgVrarIBmNOdHUFgQNk8kKovsrDrN1GIbs6kGv6sTcg4kKpH9h9HLZ72BugBQrDvCmCcRkRRtlWeYWWoimbYVjv5yoKoReWoXeu20ragLAAvCJieEDM3lvcyna1jMVPpdnZBz96mjoiIYK+oUSGdyfwSEBWGTU7Xxmyxlyklp3O+qhJyDdHSjk5Nwb4iiD8QhttEsasZsWGV0jKVMuZc3L4P9rfdJJIrsP45e/KG7tLm8iuL8XQRrIC3XivE3YkpAzEkCGW+DSmuF+cPbHmnfUfT6edsstEhZS9jkV8b74Eei2N98PzhQUY207YLbvfDUYYru97A8OyKkiuy3YUqL2xIjk7LltL5yd+UgXl68e9RZ7ghj8sULxNYt4TMHWjohOQVzU9jvonhYd4QpAAeicbc1LN594wMHWg+iAz1IRzdEimQGS8jzCdunwRSIA+JgbzmfRaJI6z4Yvo61houDtYRNAbkBOK1RImIj9f5wk8FT0wKTd92vvv9htI0+jSmoiKTcecTGJ+6MCaIxWJjFWsT+spaw8UXUWWAhXUUiMu13VULKgW31bra1eLnXVeSTkE9RsyBsfCFkiDpJ5tJ1VEQm/K5OSHmBODnt+2o6WzFnjA8ikiLCEjOpJr+rEmKOG4CtW8I31hI2vkpEpplfqmJ6qYWq2F2/qxNSjtslsTjnDtoVum0W8u4Iawkb31XHhkima2x3Dj+JA04UUotYi7iwrCVsikJj2VUGZp8j4ixaH7EvvAU1i3OQdgq8kCNYO2HkmrWETVEQMrRVnGJgeo8lhfeNN4d4MWk7chSQBWFTNCKSYmfVZa5PdVlSeN94A3XJaQrWLaG20acxRSMRmaZ92wDXJo+QVust84cXiOemCveStr2RMcWjKnaXhvIpbkw/i9pb1B8SdWdMZNJ+16Tk2TvcFKXGsqvEI8rg7HG/qxJe0TiI5P1l1LojjClO7ZWnmFmEu8luv6sSUo6/y5lDwn7CpmgJGTprznF7Os7Ywh6/q2PyxpuilusjICwIm6IWlUW66q9z/V6VLeYwJcmCsCl65ZEpDtSP0DveQjJd7Xd1TK4ptsecMcWuOjbErpp5esZ2Wy7iEqRpzfkRFBaETWA0lV2mZVuUy/e6bA6xKRkWhE2gtFeeIuo4fHDf5hCXDOuOMCY4lmdM3F9I2RxiUxIsCJvAicoi3fXXGLi/xHCyy+/qmCemtmzZmKBJRKY53DjEtXsx7i0+7Xd1jNkyG90wgbUtOkJXfSUXx2o42pRkW3TE7yqZrVDbY86YwKpL9LO7Bs6PttocYhNIFoRN4LVWnKO5Mm5ziIPMli0bE2w7q05TEYvQO2EJ4QPHsqgZE3wRSbG/9jypjHJt8ojNITaBYe9UUzKWp67dm1/ixvRzflfHbIItWzamRCxPXRu4P8/tuWN+V8eYddkUNVNytkVH6G6o5PxoFVHporm81+8qmTWoBqsPN9esJWxKUl2in4P1KS6NR5lY2OV3dYxZlbWETclqLu9lsfY5esaEY00tVMXu+l0ls4pMgPpwc81awqaktVeeorkyzrnRZlvMUaxsipoxpUvIsLf6DFXxqC3mMEXJgrApeRFJ0VXnDs5ZQvjio4BmMjk/gsKCsAmFmJOku+EDphfT3Jh+1u/qGPOABWETGuWRKQ43DnFnZoGhuSN+V8cs09wv1Ci5xRoiUiMi/0lELovIJRH5qIjUicg7InLN+1qb9fzXRaRPRK6IyKeyyo+LyAXvsa+IiOTjooxZjTuHeIK+ezC52OF3dUyREpFf8eLdeRH5jojU5Ou1NtoS/jXg+6p6ADgCXAJeA95V1U7gXe8+ItIFvAwcAl4EvioiEe/7fA14Fej0jhdzdB3GbFhNfICD9Ul6xupsxkSRKMLZEe8A3ap6GLgKvP7EF7mKdYOwiGwHfgz4TQBVXVTVSeAl4C3vaW8Bn/ZuvwR8S1UXVLUf6AOeF5EdwHZV/YGqKvCNrHOMKajGsqvsroHe8T2Wdc1vWny5I1T1j1U15d39K6D9ia9zFRtpCe8GRoH/R0TOiMhviEgl0KyqdwC8r03e89uAgazzB72yNu/2o+WPEZFXReSkiJwcHZ3a1AUZs1GtFedorIxZ1rXS1bAcR7zj1S1+n38AfC+XFcu2kbk6UeBZ4B+r6g9F5Nfwuh5WsVI/r65R/nih6pvAmwDHjj8TnB52Ezgdle9zdek4g7PH6ah83+/qhFaeFleMqeqJ1R4UkT8BWlZ46Euq+l3vOV8CUsA381FB2FgQHgQGVfWH3v3/hBuEh0Vkh6re8boaRrKenz3i0Q4MeeXtK5SvaW4pyrWp52nfdoXyiLWKTW4tL+a4fO8oo5F9NJZd9btKpkBU9RNrPS4irwB/A/i414WaF+v+D6aqd4EBEdnvFX0c6AXeBl7xyl4Bvuvdfht4WUQSIrILdwDuPa/LYlpEXvBmRXw265xVRRxh4P4879/dyeDscVvxZHJuOSH80EwV00srNYxMPqlCJqM5P56EiLwI/CLw06o6l5MLXcVGlw79Y+CbIhIHPgD+Pm4A/7aIfA64BXwGQFUvisi3cQN1CviCqqa97/N54LeActw+lnX7Wcqjc+yvj3BtIs3ViSVuRfezp2aRhrKrRCS13unGbEhUFjlYd5lLEwc4UDtLIjLtd5WMv34dSADveDNp/0pV/2E+XmhDQVhVzwIr9a18fJXnvwG8sUL5SaB7E/VDUNoqzrA91sLVe61MLaS4OBalOnGYXdX3qE3cRAjOEkVTvOLOLPtq++ib2s/+2vNEZdHvKoVE8S2uUNW9hXqtwAwJV8XucqSxh866BI7A1EKKsyNV9IwfY2qp1e/qmRJRHpni6e23uDl91GZMFIplUQuOqCzSUfk+x5tHqE64jfjRuSXODDdw6d5zzKXrfK6hKQXboiM0lA/a0mZTEIEKwsuWW8X76mI4AhmFOzMLnLzbwd1kt2XJMk+sOjZEeWSSsYU9flclFKwlHEBRWaS98hQnWu5SW+YG3VRG6R1z6J04wmyqwecamqCrS/SjGrEZEyavAhuEl22LjnC44TxdDRniEXc9yOjcEqeG2xiaO2KtYvNEGsuukkzXWI6JPNIiXLZcSIEPwuDO82wp7+F48w12bEsAbqv48rjSM36EmVTTOt/BmNU1lV1mZqnZ5qibvCiJILysPDLFgdpTHG6apTLmXtp4conTwzsYnD1uiVrMljWU9TG11GozJvJCQ72zRsn9ry5kaEhcp7qpnLvJbvonF0lllKsTS9yZOURn7Rg18YH1v5ExWYQM9Yl+phbb7P2Ta153RFiV7J/1mJP0prPdprEiBsD0YpqzI7Vcm3qehXSVzzU0QSNkqIrdse4tk1Ml1xJ+VGV0jO76Cca37eL6vSpmlzIM3J9ndG6PLX82mxaRFOWRCRbSVba0OYeCNKUs10q2JZxtuYvi2aYr7KuLEXWE+VSGi2NResaP2BQksykRSRF1kjbGYHKi5FvC2WJOkvbKU9SXVTM4s5/b0/OMJ5e4N9/E09VP01bZQ9yZ9buaJgAikkJxUBzLXfKElrOohVUoWsKPKo9M0Vn9Hseax6gti5JR6J9McnK4k+Fkl42Amw0RMhaAzRMLdbSpjg3RXX+Jqri7D+lyF8Xle5a32Gyc/dF+crZYI8RiTpLuhg8oi374o7gzs8DZ0QM2Cm42RMiwmKn0uxrBpbnPGxGkgb7QB2FwuyeeabhL1PlwG7zpxTSnh3cwOr/Px5qZoIg5SVvabLbEgrCnKnaX7ob7ZMVhUhnlwmgZ/dMfsX85zZqEDI61iLfMuiMM4GbNOlj/+Jzh/skkfVMnLBCbNSUi08yl6ixplNkUiyqPaC7vZV9d7LHygfvzFojNuqrjt5lY2O13NYLFdtYwj2qrPMPT1WWPlVsgNusRMlTHb1syeLNhFk1WIGTYWXX6QVrMbBaIzXriziwRSdlKzA1SrCVsVhCRFJ0156gvX71rwvr+zGpq4zeZWmyzgbqNsKTuZjVRWaSrrvfBpqLZBu7Pc3XymAVis6qWigsMzByy94hZkwXhdcScJN311x4kic92Z2aB3okjlsjFrCgqi7RUfMDtueN+V6XIKZlM7o+gsCC8AYnINM803npoVd2y0bklLo4/Y/92mhVVRseIO7MMJ7v8roopUhaEN6giMsHhxqEHm4lmG08ucWFsvyWKNytqLu9lYr6SqaVWv6tSlBTIZHJ/BIUF4U3YFh3hmcbRh5Y3L5taSHF2dA+zqQYfamaKmZBhT3Uvffda7A+1eYwF4U2qjg3xTOPkioF4dinDmZE2a/GYx8SdWfbW3qV3otPGEB6l1hI2m1Qbv0lXwwwrxGEW08q5kUabrG8eUx0borFCuD51xOaZP8KCsNm0hsR1DjXMrxiIUxmlZ7SSobkjha+YKWptlWeYT2XsvWEesCD8BBrLrtLVsLhiIM4oXB5X+qc/YvNEzQNChgO1V7k5pUws7PK7OkVBcT8vuT6CwoLwE2oqu8yB+tX/9+mfTHL53lHbqcM8kIhMc7D+Hr3j220g11gQzoWW8h66GlYPxMOzi5wfO2hJv80DtfGbPFVdRs9Yu80xt4E5kwvrBeKphRRnRnYxudhRwFqZYtZWcYp4xOHSxAHrsgoxC8I51FLew/76yKqPz6cynB2ptdVTBnCTRHXVXWN6McUH95/1uzq+scUaJqfaKs6smBR+WUbh4liUgdnnClgrU6wSkWkONUxye3o+vDMmrDvC5Fp75ak1AzHAtYkFm0tsALd/eE9tgqsTagt9QsiCcJ60V57iQP0Kc9ey9I5tYy5dV6AamWLWXnmK2rIYF8eaQjlQZy1hkxetFefWHKxLZZRr93bZ6imDkOFg3WUArtw7YO+JELHfdJ6tN2tiPLnE2PzeAtbIFKu4M0tX/TjjySUGZ8OTg9gG5kzetZT3cKghteLKOoArE+W2mMMAUBMfoLMuwvV7C7aiLiQsCBdIc3nvqkucF9PKjelnCl8pU5RaK87RXJmgd3x7OBb42OwIUyhNZZfpbpxdMRAP3J+3hRwGcPuH91b3EHOEyxN7Sn4hh3VHmIJqSFzncOP0ivmIr0zUW65ZA3h7GzYMMr2Y5tZMePqHw8iCsA/qEv0cbpx4LBDPLmUYmDnmU61MsamMjnGwPsnNqSSj8/v8rk7+WHeE8UNNfIBjTcOP7VnXP5m0CfvmgcayqzxdXc6l8XLLuFaiLAj7qCp2lwP1M4+VXx5vsm4J88BT205RnYjSO95Rsu8LVc35ERQWhH3WkLhOY8XDS5xnlzL03z/qT4VM0YlIigO1V1nKqG2NVILst1kE9tT0PzZjYuD+vM0TNQ+4iX5GuDOzwJ250prOaLMjjO8qIhPsrX084U/v+PZQ5hEwK6uODXGgPlN6iX5sYM4Ugx0V56iKP5yLeDGtlkfAPKSlvIeO7WVcHGtiIV3ld3VMDtinu0hEJMWBujuPdUuMzi2V3L+f5snsrDpNedShd6KzZBZyWEvYFIWq2F121z6eQ+LqhDK91OJDjUwxWt6RI5nK2EKOErBuEBaR/SJyNuu4LyI/JyJ1IvKOiFzzvtZmnfO6iPSJyBUR+VRW+XERueA99hURWTvhbgi1VZyiMvbwryWj0DveYkl+zAPLA3UD9+cZmT/gd3WeiA3MrUNVr6jqUVU9ChwH5oDvAK8B76pqJ/Cudx8R6QJeBg4BLwJfFZHlzs6vAa8Cnd7xYk6vpgREJEVX/d3HuiVmlzJcnTxk/cPmgerYEPvrlrg8nmAm1eR3dcwWbfYT/XHguqreBF4C3vLK3wI+7d1+CfiWqi6oaj/QBzwvIjuA7ar6A3VnUn8j6xyTpSp2lz21icfKh2cXw7sPmVlRc3kvHdvLOD/aGtyZNDY7YlNeBn7Hu92sqncAvK/Lf4rbgIGscwa9sjbv9qPljxGRV0XkpIicHB2d2mQVS0NrxRmqE48PulybSFu2NfOQp7adoioe4dLEgZIZqAuTDQdhEYkDPw38x/WeukKZrlH+eKHqm6p6QlVPNDaGIJ/qCiKS4mD9jceS/GQUesbqbHqSeSAiKfbXXmY+FcyMa9YnvHE/CZxW1WHv/rDXxYD3dcQrHwSym2rtwJBX3r5CuVlFRWSC/XVLj5UvppXB2YM+1MgUq7gzS3fDILenA7jS0rojNuxn+LArAuBt4BXv9ivAd7PKXxaRhIjswh2Ae8/rspgWkRe8WRGfzTrHrKK5vJcd2x7vH749vVCyyVzM1lRGx+iqv8/liWr7TylANhSERaQC+CTwe1nFXwY+KSLXvMe+DKCqF4FvA73A94EvqGraO+fzwG/gDtZdB76Xg2soeXurex6btpbKKPcWdvpTIVO06hL9PF0tXJ3cF6iZNBnN/REUG+rFV9U5oP6RsnHc2RIrPf8N4I0Vyk8C3ZuvZrjFnCSHGoY4ebfloTfXjalKGsv8q5cpTq0V55hbOsGduWdorTjnd3XMOoLzpzLktkVHOFifeqhsejFt80PNY4QMu7efZjwZDcT7wwbmTGAszwnNNjD9tE+1McVsecZE/1TpJoIvFRaEA2bX9rMPzR8enl0I7iR9k1dxZ5bd1TeLf99Cmx1hgiQqi3TVX3+wN11GYSQZ7NwBJn8qo2PUJoYYW9jjd1VWZd0RJnDKI1M80zj6IL/EjamUrZQyq6qJD6AaYS5d53dVAkNE/hcROe8lLftjEclbFn0LwgFVHRt6MFC3mFbGF/b6XCNTzBrLrjKz1FS0f6yLsCX8K6p62Etc9vvALz3xd1yFBeEAay7vZVeNm96yf7IsUPNCTeE1ll1lyvKObIiq3s+6W8kqKRZyQYp9a2gRmQau+F0PHzQAY35XwidhvfZivO6nVbUxny8gIt/HvfZcKwPms+6/qapvbvRkEXkDd2XvFPDXVHU0x/VzXycAQfikqp7wux6FFtbrhvBee1iv2y8i8ifASlvWfElVv5v1vNeBMlX95XzUozg7iIwxJs9U9RMbfOpvA38A5CUIWyeiMcY8QkQ6s+7+NHA5X68VhJbwhvtwSkxYrxvCe+1hve5i9GUR2Q9kgJvAP8zXCxV9n7AxxpQy644wxhgfWRA2xhgfFW0QFpEXReSKiPSJyGt+1ycXROTrIjIiIj1ZZXUi8o6IXPO+1mY99rp3/VdE5FNZ5cdF5IL32Fe8nUqKloh0iMificglEbkoIl/0ykv62kWkTETeE5Fz3nX/S6+8pK/bbJKqFt0BRHB33tgNxIFzQJff9crBdf0Y8CzQk1X2vwOvebdfA/6Vd7vLu+4EsMv7eUS8x94DPoq7eer3gJ/0+9rWue4dwLPe7Srgqnd9JX3tXh23ebdjwA+BF0r9uu3Y3FGsLeHngT5V/UBVF4FvAS/5XKcnpqp/AUw8UvwS8JZ3+y3g01nl31LVBVXtx90S6nlvU9XtqvoDdT+d38g6pyip6h1VPe3dngYuAW2U+LWra8a7G/MOpcSv22xOsQbhNmAg6/6gV1aKmtXdBBXv6/JWCKv9DNq824+WB4KI7ASO4bYKS/7aRSQiImdxdyN/R1VDcd1m44o1CK/U3xW2uXSr/QwC+7MRkW3A7wI/pw8nSHnsqSuUBfLaVTWtbiaudtxW7Vp7LJbMdZuNK9YgPAhkp3tqB4Z8qku+DXv/buJ9HfHKV/sZDHq3Hy0vaiISww3A31TV5V27Q3HtAKo6Cfw58CIhum6zvmINwu8DnSKyS0TiwMvA2z7XKV/eBl7xbr8CfDer/GURSYjILqATeM/793VaRF7wRsg/m3VOUfLq+ZvAJVX91ayHSvraRaRRRGq82+XAJ3CXv5b0dZtN8ntkcLUD+CncUfTruFmNfK9TDq7pd4A7wBJu6+ZzQD3wLnDN+1qX9fwvedd/hazRcOAE0OM99ut4Kx+L9QB+FPff5/PAWe/4qVK/duAwcMa77h7gl7zykr5uOzZ32LJlY4zxUbF2RxhjTChYEDbGGB9ZEDbGGB9ZEDbGGB9ZEDbGGB9ZEDbGGB9ZEDbGGB/9//ZAqrFJAcwwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "vm = 3\n",
    "fig, ax = plt.subplots(figsize=(7, 7))\n",
    "axim = ax.imshow(dem_sardem - dem_isce, vmin=-vm, vmax=vm, cmap=\"RdYlBu\")\n",
    "cbar = fig.colorbar(axim)\n",
    "cbar.set_label(\"[meters]\")\n",
    "ax.set_title(\"sardem - isce2\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "098dfb59-0f36-46ae-a4c2-c35c65cda334",
   "metadata": {},
   "source": [
    "The lines of +/-1 difference are due to integer truncation differenceds, and differences in the EGM96 -> WGS84 geoid conversion (since `sardem` uses GDAL for this step)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d08ffc3-dc6d-4157-930d-d24d4cb245b5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}