Satpy bilinear

satpy_bilinear.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 25,
   "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": 5,
   "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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S3B_SL_1_RBT____20190211T175740_20190211T180040_20190211T195003_0179_022_041_2340_LN2_O_NR_003.SEN3\n"
     ]
    }
   ],
   "source": [
    "base_dir = Path(r'D:\\Data\\Sentinel\\Sentinel3\\SLSTR')\n",
    "print(*map(lambda x: x.name, filter(lambda x: x.name.endswith('SEN3'), base_dir.glob('*'))), sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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=\"slstr\",\n",
    "                               start_time=start_time,\n",
    "                               end_time=end_time,\n",
    "                               base_dir=os.fspath(base_dir),\n",
    "                               reader='nc_slstr')\n",
    "\n",
    "print(len(files))\n",
    "scn = Scene(filenames=files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "bands = ['S3_an', 'S2_an', 'S1_an', 'longitude_an', 'latitude_an']\n",
    "scn.load(bands)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resample to area: Satpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for reduce in tqdm([True, False]):\n",
    "    \n",
    "    for resolution in tqdm([500, 750, 1000, 1500], 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",
    "        #     new_scn.save_datasets(filename=f'SLSTR_{start_time.strftime(\"%Y%m%d\")}_{resample_alg}_{xres}m_raw.png')\n",
    "        #     new_scn.save_datasets(filename=f'SLSTR_{start_time.strftime(\"%Y%m%d\")}_{resample_alg}_{xres}m_raw.tif', writer='GeoTIFF')\n",
    "\n",
    "            compositor = satpy.composites.GenericCompositor(\"quicklook\")\n",
    "            composite = compositor([new_scn['S3_an'],\n",
    "                                    new_scn['S2_an'],\n",
    "                                    new_scn['S1_an']])\n",
    "\n",
    "            img = satpy.writers.to_image(composite)\n",
    "            img.stretch(\"linear\")\n",
    "            img.gamma(2.0)\n",
    "            img.save(f'SLSTR_{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": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {