Satpy bilinear olci

satpy_bilinear_olci.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from satpy import Scene, find_files_and_readers\n",
    "import satpy\n",
    "from pyresample.geometry import AreaDefinition, SwathDefinition\n",
    "from pyresample import bilinear, kd_tree\n",
    "\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import numpy as np\n",
    "\n",
    "import datetime\n",
    "from pathlib import Path\n",
    "import os\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# from satpy.utils import debug_on; debug_on()\n",
    "# print(satpy.config.get_environ_config_dir())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'date': '2019-01-28T10:47:01-0600',\n",
       " 'dirty': False,\n",
       " 'error': None,\n",
       " 'full-revisionid': 'ee1ee7deb9f7715a7d11cb32c6bdb6339986c08a',\n",
       " 'version': '0.11.2'}"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "satpy.version.get_versions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_area(resolution):\n",
    "    # generate a target area with \"resolution\"\n",
    "    \n",
    "    proj = dict(proj='aea', lat_1='29.5', lat_2='45.5', lat_0='37.5', lon_0='-96', x_0='0', y_0='0', datum='NAD83', units='m', no_defs='')\n",
    "    ext = [-2300000, -300000, -1200000, 600000]\n",
    "    xres, yres = resolution, -resolution\n",
    "\n",
    "    xs = abs((ext[2] - ext[0]) // xres)\n",
    "    ys = abs((ext[3] - ext[1]) // yres)\n",
    "\n",
    "    name = 'AEA_cal'\n",
    "    area_slstr = AreaDefinition(name, name, name, proj, xs, ys, ext)\n",
    "    \n",
    "    return area_slstr\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load scene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S3B_OL_1_EFR____20190211T175740_20190211T180040_20190211T193826_0179_022_041_2340_LN1_O_NR_002.SEN3\n"
     ]
    }
   ],
   "source": [
    "base_dir = Path(r'D:\\Data\\Sentinel\\Sentinel3\\OLCI')\n",
    "print(*map(lambda x: x.name, filter(lambda x: x.name.endswith('SEN3'), base_dir.glob('*'))), sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n"
     ]
    }
   ],
   "source": [
    "start_time = datetime.datetime(2019, 2, 11, 0, 0)\n",
    "end_time = start_time + datetime.timedelta(days=1)\n",
    "\n",
    "files = find_files_and_readers(sensor=\"olci\",\n",
    "                               start_time=start_time,\n",
    "                               end_time=end_time,\n",
    "                               base_dir=os.fspath(base_dir),\n",
    "                               reader='olci_l1b')\n",
    "\n",
    "print(len(files))\n",
    "scn = Scene(filenames=files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bands = ['S3_an', 'S2_an', 'S1_an', 'longitude_an', 'latitude_an']\n",
    "bands = ['Oa08', 'Oa06', 'Oa04', 'longitude', 'latitude']\n",
    "scn.load(bands)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resample to area: Satpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4a6520fea3c6457fb00e2fb44f752566",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=2), 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=(IntProgress(value=0, max=3), 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=(IntProgress(value=0, max=2), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Miniconda3\\envs\\py36\\lib\\site-packages\\dask\\core.py:137: RuntimeWarning: invalid value encountered in less\n",
      "  return func(*args2)\n",
      "C:\\Miniconda3\\envs\\py36\\lib\\site-packages\\dask\\core.py:137: RuntimeWarning: invalid value encountered in greater\n",
      "  return func(*args2)\n",
      "C:\\Miniconda3\\envs\\py36\\lib\\site-packages\\pyresample\\bilinear\\xarr.py:212: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  res = data.values[slices]\n",
      "C:\\Miniconda3\\envs\\py36\\lib\\site-packages\\pyresample\\bilinear\\xarr.py:213: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  res[mask_slices] = fill_value\n",
      "C:\\Miniconda3\\envs\\py36\\lib\\site-packages\\dask\\array\\core.py:3288: RuntimeWarning: invalid value encountered in greater\n",
      "  result = function(*args, **kwargs)\n",
      "C:\\Miniconda3\\envs\\py36\\lib\\site-packages\\dask\\array\\core.py:3288: RuntimeWarning: invalid value encountered in less\n",
      "  result = function(*args, **kwargs)\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=2), 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=(IntProgress(value=0, max=2), 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=(IntProgress(value=0, max=3), 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=(IntProgress(value=0, max=2), 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=(IntProgress(value=0, max=2), 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=(IntProgress(value=0, max=2), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "for reduce in tqdm([True, False]):\n",
    "    \n",
    "    for resolution in tqdm([300, 500, 750], leave=False):\n",
    "        \n",
    "        for resample_alg in tqdm(['nearest', 'bilinear'], leave=False):\n",
    "\n",
    "            area_slstr = get_area(resolution)\n",
    "            new_scn = scn.resample(area_slstr, resampler=resample_alg, nprocs=1, radius=10000, reduce_data=reduce)\n",
    "\n",
    "            compositor = satpy.composites.GenericCompositor(\"quicklook\")\n",
    "            composite = compositor([new_scn['Oa08'],\n",
    "                                    new_scn['Oa06'],\n",
    "                                    new_scn['Oa04']])\n",
    "\n",
    "            img = satpy.writers.to_image(composite)\n",
    "            img.stretch(\"linear\")\n",
    "            img.gamma(2.0)\n",
    "\n",
    "            img.save(f'OLCI_{start_time.strftime(\"%Y%m%d\")}_reduce_{reduce}_{resolution}m_{resample_alg}.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {