analemma.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mlubej/.pyenv/versions/3.7.6/envs/base/lib/python3.7/site-packages/bolt/utils.py:3: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated since Python 3.3,and in 3.9 it will stop working\n",
      "  from collections import Iterable\n",
      "/home/mlubej/.pyenv/versions/3.7.6/envs/base/lib/python3.7/site-packages/bolt/factory.py:17: DeprecationWarning: inspect.getargspec() is deprecated since Python 3.0, use inspect.signature() or inspect.getfullargspec()\n",
      "  args = inspect.getargspec(func)\n"
     ]
    }
   ],
   "source": [
    "# auto reload\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import os\n",
    "import imageio\n",
    "import utm\n",
    "import numpy as np\n",
    "from tqdm.auto import tqdm\n",
    "from concurrent.futures import ProcessPoolExecutor\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "from eolearn.core import FeatureType, LinearWorkflow, EOExecutor, EOTask, SaveTask, LoadTask, OverwritePermission\n",
    "from eolearn.mask import AddValidDataMaskTask\n",
    "from eolearn.io import SentinelHubInputTask\n",
    "from eolearn.coregistration import ThunderRegistration\n",
    "\n",
    "from sentinelhub import DataSource, BBox, CRS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Task definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ValidData:\n",
    "    def __init__(self, is_data_mask, clm_mask):\n",
    "        self.is_data_mask = is_data_mask\n",
    "        self.clm_mask = clm_mask\n",
    "        \n",
    "    def __call__(self, eopatch):\n",
    "        return np.logical_and(eopatch.mask[self.is_data_mask].astype(np.bool),\n",
    "                              np.logical_not(eopatch.mask[self.clm_mask].astype(np.bool)))\n",
    "    \n",
    "\n",
    "class ValidCoverageTask(EOTask):\n",
    "    def __init__(self, feature_in, feature_out):\n",
    "        self.feature_in = feature_in\n",
    "        self.feature_out = feature_out\n",
    "        \n",
    "    def execute(self, eopatch):\n",
    "        mask = eopatch[self.feature_in]\n",
    "        coverage = np.count_nonzero(mask == 1, axis=(1,2))/np.prod(mask.shape[1:])\n",
    "        eopatch[self.feature_out] = coverage\n",
    "        return eopatch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download RBG L2A bands for plotting\n",
    "download_task = SentinelHubInputTask(\n",
    "    bands_feature=(FeatureType.DATA, 'RGB'),\n",
    "    bands = ['B04', 'B03', 'B02'],\n",
    "    resolution=10,\n",
    "    maxcc=1.0,\n",
    "    time_difference=timedelta(hours=2),\n",
    "    data_source=DataSource.SENTINEL2_L2A,\n",
    "    max_threads=10,\n",
    "    additional_data=[\n",
    "        (FeatureType.MASK, 'dataMask'),\n",
    "        (FeatureType.MASK, 'CLM'),\n",
    "        (FeatureType.DATA, 'CLP'),\n",
    "        (FeatureType.DATA, 'sunZenithAngles'),\n",
    "        (FeatureType.DATA, 'sunAzimuthAngles')\n",
    "    ])\n",
    "\n",
    "# task for creating a valid data mask\n",
    "valid_data_task = AddValidDataMaskTask(ValidData(is_data_mask='dataMask', clm_mask='CLM'), 'VALID_DATA')\n",
    "        \n",
    "\n",
    "# task for calculating the valid coverage\n",
    "valid_cov_task = ValidCoverageTask(feature_in=(FeatureType.MASK, 'VALID_DATA'), \n",
    "                                   feature_out=(FeatureType.SCALAR, 'VALID_COVERAGE'))\n",
    "\n",
    "\n",
    "# task for coregistrating the time frames\n",
    "coreg_task = ThunderRegistration((FeatureType.DATA, 'RGB'), valid_mask_feature = (FeatureType.MASK, 'VALID_DATA'), \n",
    "                                 channel=0)\n",
    "\n",
    "\n",
    "# tasks for loading and saving\n",
    "load_task = LoadTask('./eopatches', lazy_loading=True)\n",
    "save_task = SaveTask('./eopatches', overwrite_permission=OverwritePermission.OVERWRITE_PATCH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Download"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#(lat, lon), width in meters\n",
    "location_data = [\n",
    "    [[30.962476, 34.730068], 1000], # ashalim\n",
    "    [[25.197020, 55.274212], 1250], # burj khalifa\n",
    "    [[29.979221, 31.134213], 1e3], # pyramids\n",
    "    [[35.710054, 139.810714],1250] # tokyo sky tree\n",
    "]\n",
    "\n",
    "factors = [2.0, 2.5, 2.5, 2.75]\n",
    "th = 0.9\n",
    "\n",
    "time_interval = [datetime(2019,1,1), datetime(2019,12,31)]\n",
    "bbox_list = []\n",
    "for idx in range(len(location_data)):\n",
    "    x,y,zone,letter = utm.from_latlon(*location_data[idx][0])\n",
    "    x,y = np.round([x,y], -1)\n",
    "    d = location_data[idx][1]\n",
    "    bbox_list.append(BBox((x-d, y-d, x+d, y+d), eval(f'CRS.UTM_{zone}{\"N\" if letter >= \"N\" else \"S\"}')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ec0542b8edd74e419c0ba6078f8900c6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=4.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mlubej/.pyenv/versions/3.7.6/envs/base/lib/python3.7/site-packages/holoviews/operation/datashader.py:3: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated since Python 3.3,and in 3.9 it will stop working\n",
      "  from collections import Callable, Iterable\n"
     ]
    }
   ],
   "source": [
    "os.system('rm -rf eopatches && mkdir eopatches')\n",
    "\n",
    "workflow = LinearWorkflow(\n",
    "    download_task,\n",
    "    valid_data_task,\n",
    "    valid_cov_task,\n",
    "    coreg_task,\n",
    "    save_task\n",
    ")\n",
    "\n",
    "execution_args = []\n",
    "for idx, bbox in enumerate(bbox_list):\n",
    "    execution_args.append({\n",
    "        download_task: {'bbox': bbox, 'time_interval': time_interval},\n",
    "        save_task: {'eopatch_folder': f'eopatch_{idx}'}\n",
    "    })\n",
    "    \n",
    "executor = EOExecutor(workflow, execution_args, save_logs=True)\n",
    "executor.run(workers=8, multiprocess=False)\n",
    "\n",
    "executor.make_report()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# load the previously created eopatches\n",
    "\n",
    "workflow = LinearWorkflow(\n",
    "    load_task\n",
    ")\n",
    "\n",
    "execution_args = []\n",
    "for idx, bbox in enumerate(bbox_list):\n",
    "    execution_args.append({\n",
    "        load_task: {'eopatch_folder': f'eopatch_{idx}'}\n",
    "    })\n",
    "    \n",
    "eopatches = []\n",
    "for args in execution_args:\n",
    "    eopatches.append(workflow.execute(args).eopatch())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b64bf8e9ce564a4bb09dca5fa3775563",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=4.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=71.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=73.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=72.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=72.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# create RGB animations \n",
    "os.system('rm -rf graphs && mkdir graphs')\n",
    "\n",
    "def plot_image(idx, f):\n",
    "    fig = plt.figure(figsize=(10,10))\n",
    "    plt.imshow(np.clip(eop.data['RGB'][idx]*f,0,1))\n",
    "    plt.xticks([])\n",
    "    plt.yticks([])\n",
    "    plt.axis('off')\n",
    "    plt.savefig(f'graphs/true_color_{name}_{idx}.png', dpi=50, bbox_inches='tight')\n",
    "    plt.close()\n",
    "    \n",
    "for idx in tqdm(range(len(eopatches)), total=len(eopatches)):\n",
    "    eop = eopatches[idx]\n",
    "    f = factors[idx]\n",
    "    name = f'loc{idx}'\n",
    "    \n",
    "    def plot(idx):\n",
    "        plot_image(idx,f)\n",
    "\n",
    "    with ProcessPoolExecutor(max_workers=8) as executor:\n",
    "        _ = list(tqdm(executor.map(plot, range(len(eop.timestamp))), total=len(eop.timestamp), leave=False))\n",
    "\n",
    "    n_valid = np.count_nonzero(eop.scalar['VALID_COVERAGE']  > th)\n",
    "    with imageio.get_writer(f'figs/true_color_{name}.gif', mode='I', duration=2/n_valid) as writer:\n",
    "        for i in range(len(eop.timestamp)):\n",
    "            if eop.scalar['VALID_COVERAGE'][i] > th:\n",
    "                image = imageio.imread(f'graphs/true_color_{name}_{i}.png')\n",
    "                writer.append_data(image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create graph animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define function for plotting time frame\n",
    "\n",
    "def plot_graph(idx, f, name):\n",
    "    x = np.mean(eop.data['sunAzimuthAngles'][...,0], axis=(1,2))\n",
    "    y = np.mean(eop.data['sunZenithAngles'][...,0], axis=(1,2))\n",
    "    valid_mask = eop.scalar['VALID_COVERAGE'].squeeze() > th\n",
    "    ids = np.array(range(len(valid_mask)))\n",
    "    dates = np.array([ts.date().isoformat()[:-3] for ts in eop.timestamp])\n",
    "    transition_dates = [datetime(2019,3,20), datetime(2019,6,21),datetime(2019,9,23), datetime(2019,12,22)]\n",
    "    transition_doys = [(x-datetime(2019,1,1)).days for x in transition_dates]    \n",
    "    try:\n",
    "        last_valid = ids[:idx+1][valid_mask[:idx+1]][-1]\n",
    "    except:\n",
    "        last_valid = ids[valid_mask][-1]\n",
    "    \n",
    "    fig, axs = plt.subplots(2,2,figsize=(10,10))\n",
    "    plot_ids = np.array(range(len(valid_mask[:idx])))[valid_mask[:idx]].astype(int)\n",
    "\n",
    "    ax1 = axs[0,0]\n",
    "    ax1.plot(x,y)\n",
    "    ax1.plot(x[:idx+1],y[:idx+1],'r')\n",
    "    ax1.plot(x[plot_ids],y[plot_ids],'kx')\n",
    "    ax1.plot(x[idx], y[idx], 'r', marker='o')\n",
    "    ax1.set_xlabel('Azimuth [°]')\n",
    "    ax1.set_ylabel('Zenith [°]')\n",
    "\n",
    "    def idx_to_days(x):\n",
    "        return x*365/len(dates)\n",
    "\n",
    "    def days_to_idx(x):\n",
    "        return x*len(dates)/365\n",
    "    \n",
    "    ax2 = axs[0,1]\n",
    "    ax2.plot(y)\n",
    "    ax2.plot(range(len(y[:idx+1])), y[:idx+1],'r')\n",
    "    ax2.plot(np.array(range(len(y)))[plot_ids], y[plot_ids],'kx')\n",
    "    ax2.plot(idx, y[idx], 'r', marker='o')\n",
    "    ax2.set_xlabel('Time')\n",
    "    ax2.set_xticks([days_to_idx(doy) for doy in transition_doys])\n",
    "    ax2.set_xticklabels([dt.strftime(format='%Y-%m') for dt in transition_dates], ha='right')\n",
    "    ax2.set_yticks([])\n",
    "    \n",
    "    ax3 = axs[1,0]\n",
    "    ax3.plot(x, range(len(x)))\n",
    "    ax3.plot(x[:idx+1], range(len(x[:idx+1])), 'r')\n",
    "    ax3.plot(x[plot_ids], np.array(range(len(x)))[plot_ids], 'kx')\n",
    "    ax3.plot(x[idx], idx, 'r', marker='o')\n",
    "    ax3.set_xticks([])\n",
    "    ax3.set_yticks([])\n",
    "    \n",
    "    sax3 = ax3.secondary_yaxis('right', functions=(idx_to_days, days_to_idx))\n",
    "    sax3.set_ylabel('Time')\n",
    "    sax3.set_yticks(transition_doys)\n",
    "    sax3.set_yticklabels([dt.strftime(format='%Y-%m') for dt in transition_dates], rotation=90, va='top')\n",
    "\n",
    "    ax1.axvline(x=x[idx],ymin=-1.2,ymax=1,c=\"gray\",linewidth=1, linestyle='dashed',zorder=0, clip_on=False)\n",
    "    ax3.axvline(x=x[idx],ymin=0,ymax=1,c=\"gray\",linewidth=1, linestyle='dashed',zorder=0, clip_on=False)\n",
    "    ax1.axhline(y=y[idx],xmin=0,xmax=1.2,c=\"gray\",linewidth=1, linestyle='dashed',zorder=0, clip_on=False)\n",
    "    ax2.axhline(y=y[idx],xmin=0,xmax=1,c=\"gray\",linewidth=1, linestyle='dashed',zorder=0, clip_on=False)\n",
    "    \n",
    "    ax4 = axs[1,1]\n",
    "    \n",
    "    if last_valid is not None:\n",
    "        ax4.imshow(np.clip(eop.data['RGB'][last_valid]*f,0,1))\n",
    "    else:\n",
    "        pass\n",
    "    \n",
    "    ax4.set_xticks([])\n",
    "    ax4.set_yticks([])\n",
    "    ax4.axis('off')  \n",
    "        \n",
    "    # comment out below to plot\n",
    "    plt.savefig(f'graphs/graph_{name}_{idx}.png', dpi=100, bbox_inches='tight')\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f5dd24fb4b5e41aa821a22525d6f30d5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=4.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=71.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=73.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=72.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=72.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# create RGB and solar angle animation\n",
    "for idx in tqdm(range(len(eopatches)), total=len(eopatches)):\n",
    "    eop = eopatches[idx]\n",
    "    name = f'loc{idx}'\n",
    "    f = factors[idx]\n",
    "\n",
    "    def plot(idx):\n",
    "        plot_graph(idx, f, name)\n",
    "\n",
    "    os.system('rm -rf graphs && mkdir graphs')\n",
    "    with ProcessPoolExecutor(max_workers=8) as executor:\n",
    "        _ = list(tqdm(executor.map(plot, range(len(eop.timestamp))), total=len(eop.timestamp), leave=False))\n",
    "\n",
    "    with imageio.get_writer(f'figs/graph_{name}.gif', mode='I', duration=4.0/len(eop.timestamp)) as writer:\n",
    "        for i in range(len(eop.timestamp)):\n",
    "            image = imageio.imread(f'graphs/graph_{name}_{i}.png')\n",
    "            writer.append_data(image)"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}