data_sources_explorer.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data-sources exploration using `eo-learn`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook shows some examples on how to retrieve EO and non-EO data using `eo-learn`. \n",
    "\n",
    "The steps are as follow:\n",
    " * split area of interest into easy-to-process EOPatches\n",
    " * add Sentinel-2 imaging data\n",
    " * add Sentinel-1 imaging data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add generic packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "from pathlib import Path\n",
    "\n",
    "from matplotlib import dates\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import matplotlib.pyplot as plt\n",
    "import geopandas as gpd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set path to data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_dir = Path('data')\n",
    "os.listdir(data_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`eo-learn` and `sentinelhub` imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from eolearn.core import FeatureType\n",
    "from eolearn.io import SentinelHubInputTask, SentinelHubEvalscriptTask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sentinelhub import BBoxSplitter, CRS, DataCollection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Split country into smaller bounding boxes <a id='splitter'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load shapefile of Denmark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "country_filename = data_dir / 'svn_border.geojson'\n",
    "country = gpd.read_file(country_filename)\n",
    "\n",
    "country.plot()\n",
    "country.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set CRS to UTM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "country_crs = CRS.UTM_32N\n",
    "country = country.to_crs(country_crs.pyproj_crs())\n",
    "\n",
    "country.plot()\n",
    "country.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get size of country in pixels to decide number and size of bounding boxes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "country_shape = country.geometry.values[-1]\n",
    "\n",
    "width_pix = int((country_shape.bounds[2] - country_shape.bounds[0]) / 10)\n",
    "height_pix = int((country_shape.bounds[3] - country_shape.bounds[1]) / 10)\n",
    "\n",
    "print(f'Dimension of the area is {width_pix} x {height_pix} pixels')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Split area into 45x35 boxes bounding "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "bbox_splitter = BBoxSplitter([country_shape], country_crs, (45, 35))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "geometry = [bbox.geometry for bbox in bbox_splitter.get_bbox_list()]\n",
    "bbox_list = bbox_splitter.get_bbox_list()\n",
    "idxs_x = [info['index_x'] for info in bbox_splitter.get_info_list()]\n",
    "idxs_y = [info['index_y'] for info in bbox_splitter.get_info_list()]\n",
    "\n",
    "gdf = gpd.GeoDataFrame(\n",
    "    {'index_x':idxs_x, 'index_y':idxs_y},\n",
    "    geometry=[bbox.geometry for bbox in bbox_list],\n",
    "    crs=bbox_list[0].crs.pyproj_crs()\n",
    ")\n",
    "\n",
    "gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# if bboxes have all same size, estimate offset\n",
    "xl, yl, xu, yu = gdf.geometry[0].bounds\n",
    "xoff, yoff = (xu - xl) / 3, (yu - yl) / 5\n",
    "\n",
    "# figure\n",
    "fig, ax = plt.subplots(figsize=(45,35))\n",
    "gdf.plot(ax=ax, facecolor='w', edgecolor='r', alpha=0.5, linewidth=2)\n",
    "country.plot(ax=ax, facecolor='w', edgecolor='b', alpha=0.5, linewidth=2.5)\n",
    "ax.set_title('Denmark tiled in a 45 x 35 grid');\n",
    "\n",
    "# add annotiation text\n",
    "fontdict = {'family': 'monospace', 'weight': 'normal', 'size': 14}\n",
    "for idx in gdf.index:\n",
    "    eop_name = '{0}x{1}'.format(gdf.index_x[idx], gdf.index_y[idx])\n",
    "    centroid, = list(gdf.geometry[idx].centroid.coords)\n",
    "    ax.text(centroid[0] - xoff, centroid[1] + yoff, str(idx), fontdict=fontdict)\n",
    "    ax.text(centroid[0] - xoff, centroid[1] - yoff, eop_name, fontdict=fontdict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Retrieve S2 L1C data <a id=\"sentinel-2\"></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "s2_rgb_task = SentinelHubInputTask(\n",
    "    data_collection=DataCollection.SENTINEL2_L1C,\n",
    "    bands=['B04', 'B03', 'B02'],\n",
    "    bands_feature=(FeatureType.DATA, 'S2-RGB'),\n",
    "    additional_data=[(FeatureType.MASK, 'dataMask')],\n",
    "    resolution=(10, 10),\n",
    "    maxcc=0.1\n",
    ")\n",
    "\n",
    "ndvi_evalscript = \"\"\"\n",
    "//VERSION=3\n",
    "\n",
    "function setup() {\n",
    "  return {\n",
    "    input: [\"B04\", \"B08\"],\n",
    "    output:[\n",
    "      {\n",
    "        id: \"ndvi\",\n",
    "        bands: 1,\n",
    "        sampleType: SampleType.FLOAT32\n",
    "      },\n",
    "    ]\n",
    "  }\n",
    "}\n",
    "\n",
    "function evaluatePixel(sample) {\n",
    "  let ndvi = index(sample.B08, sample.B04);\n",
    "  return {\n",
    "    ndvi: [ndvi],\n",
    "  };\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "s2_ndvi_task = SentinelHubEvalscriptTask(\n",
    "    features=[(FeatureType.DATA, 'ndvi', 'NDVI')],\n",
    "    evalscript=ndvi_evalscript,\n",
    "    data_collection=DataCollection.SENTINEL2_L1C,\n",
    "    resolution=(10, 10),\n",
    "    maxcc=0.1\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "time_interval = ['2019-05-01','2019-09-01']\n",
    "idx = 436\n",
    "bbox = bbox_splitter.bbox_list[idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download TRUE-COLOR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "eopatch = s2_rgb_task.execute(bbox=bbox, time_interval=time_interval)\n",
    "\n",
    "eopatch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "eopatch.timestamp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download NDVI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "eopatch = s2_ndvi_task.execute(eopatch)\n",
    "\n",
    "eopatch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "eopatch.timestamp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot RGB of time frames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "time_idx = 3\n",
    "\n",
    "rgb = eopatch.data['S2-RGB']\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(15, 15))\n",
    "ax.imshow(rgb[time_idx] * 3.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the median RGB values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15, 15))\n",
    "\n",
    "ax.imshow(np.median(rgb, axis=0).squeeze() * 3.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the median NDVI values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "eopatch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ndvi = eopatch.data['NDVI']\n",
    "median_ndvi = np.median(ndvi, axis=0).squeeze()\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(15,15))\n",
    "im = ax.imshow(median_ndvi, cmap=plt.cm.YlGn)\n",
    "\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes('right', size='5%', pad=0.05)\n",
    "fig.colorbar(im, cax=cax, orientation='vertical');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot temporal NDVI of a given location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dates_num = dates.date2num(eopatch.timestamp)\n",
    "dates_str = [timestamp.date().isoformat() for timestamp in eopatch.timestamp]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(15, 15))\n",
    "ax.plot(dates_num, ndvi[:, 100, 520, :].squeeze(), 'g')\n",
    "\n",
    "ax.set_title('NDVI evolution')\n",
    "ax.set_xticks(dates_num);\n",
    "ax.set_xticklabels(dates_str, rotation=45, ha='right');\n",
    "ax.set_ylabel('NDVI');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Retrieve S1 data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "s1_iw_des_task = SentinelHubInputTask(\n",
    "    data_collection=DataCollection.SENTINEL1_IW_DES,\n",
    "    bands=['VV'],\n",
    "    bands_feature=(FeatureType.DATA, 'S1-IW-DES'),\n",
    "    additional_data=[(FeatureType.MASK, 'dataMask')],\n",
    "    resolution=(10, 10)\n",
    ")\n",
    "\n",
    "s1_iw_asc_task = SentinelHubInputTask(\n",
    "    data_collection=DataCollection.SENTINEL1_IW_ASC,\n",
    "    bands=['VV'],\n",
    "    bands_feature=(FeatureType.DATA, 'S1-IW-ASC'),\n",
    "    additional_data=[(FeatureType.MASK, 'dataMask')],\n",
    "    resolution=(10, 10)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "eopatch_s1_des = s1_iw_des_task.execute(bbox=bbox, time_interval=['2019-07-01','2019-08-01'])\n",
    "\n",
    "eopatch_s1_des"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vv_data = eopatch_s1_des.data['S1-IW-DES']\n",
    "vv_data[np.isnan(vv_data)] = 0\n",
    "\n",
    "vv_des_r = np.percentile(vv_data, 80, axis=0)\n",
    "vv_des_g = np.percentile(vv_data, 20, axis=0)\n",
    "vv_des_b = np.std(vv_data, axis=0)\n",
    "\n",
    "vv_rgb = np.concatenate((vv_des_r, vv_des_g, vv_des_b), axis=-1)\n",
    "\n",
    "plt.figure(figsize=(15, 15))\n",
    "plt.imshow(vv_rgb);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "eopatch_s1_asc = s1_iw_asc_task.execute(bbox=bbox, time_interval=['2019-07-01','2019-08-01'])\n",
    "\n",
    "eopatch_s1_asc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vv_data = eopatch_s1_asc.data['S1-IW-ASC']\n",
    "vv_data[np.isnan(vv_data)] = 0\n",
    "\n",
    "vv_des_r = vv_data[0]\n",
    "vv_des_g = vv_data[2]\n",
    "vv_des_b = vv_data[10]\n",
    "\n",
    "vv_rgb = np.concatenate((vv_des_r, vv_des_g, vv_des_b), axis=-1)\n",
    "\n",
    "plt.figure(figsize=(15, 15))\n",
    "plt.imshow(vv_rgb);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, Sentinel-2 L2A data can be added, as well as Digital Elevation data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "eo",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}