use STAC metadata to create xr object

STAC_xr.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "inner-prediction",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import satsearch\n",
    "from affine import Affine\n",
    "import numpy as np\n",
    "from dask import array as darray\n",
    "import dask\n",
    "import datetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "single-performance",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "500"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "URL='https://earth-search.aws.element84.com/v0'\n",
    "results = satsearch.Search.search(url=URL,\n",
    "                                  collections=['sentinel-s2-l2a-cogs'],\n",
    "                                  bbox=[-104.791,39.783,-104.677,40.170],\n",
    "                                  datetime='2017-06-27/2021-02-27')\n",
    "stac_items = results.items().geojson()['features']\n",
    "len(stac_items)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "occupational-fiber",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _delayed_open(stac_item,band):\n",
    "    transform = stac_item['assets'][band]['proj:transform'][0:6]\n",
    "    nx,ny = stac_item['assets'][band]['proj:shape']\n",
    "    x, y = Affine(*transform) * (np.arange(nx)+0.5, np.arange(ny)+0.5)\n",
    "    href = stac_item['assets'][band]['href']\n",
    "    date = datetime.datetime.fromisoformat(stac_item['properties']['datetime'][0:10]).date()\n",
    "    arr = darray.from_delayed(dask.delayed(xr.open_rasterio)(href).data, (1,nx,ny), dtype=float)[np.newaxis,:,:,:]\n",
    "    return(xr.DataArray(data=arr, coords=([date],[band],y,x),dims=('time','band','y','x')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "compliant-connecticut",
   "metadata": {},
   "outputs": [],
   "source": [
    "def stack_stac(stac_items,asset):\n",
    "    return(xr.concat([_delayed_open(stac_item,'B01') for stac_item in stac_items],dim='time').sortby('time'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "popular-diversity",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "836 ms ± 44.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "stack_stac(stac_items,'B01')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "norwegian-password",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>&lt;xarray.DataArray &#x27;getitem-fac2a9b2409c6bece12610e9e599b4a3&#x27; (time: 500, band: 1, y: 1830, x: 1830)&gt;\n",
       "dask.array&lt;getitem, shape=(500, 1, 1830, 1830), dtype=float64, chunksize=(1, 1, 1830, 1830), chunktype=numpy.ndarray&gt;\n",
       "Coordinates:\n",
       "  * time     (time) object 2017-06-28 2017-07-01 ... 2021-02-25 2021-02-27\n",
       "  * band     (band) &lt;U3 &#x27;B01&#x27;\n",
       "  * y        (y) float64 4.5e+06 4.5e+06 4.5e+06 ... 4.39e+06 4.39e+06 4.39e+06\n",
       "  * x        (x) float64 5e+05 5.001e+05 5.001e+05 ... 6.097e+05 6.098e+05</pre>"
      ],
      "text/plain": [
       "<xarray.DataArray 'getitem-fac2a9b2409c6bece12610e9e599b4a3' (time: 500, band: 1, y: 1830, x: 1830)>\n",
       "dask.array<getitem, shape=(500, 1, 1830, 1830), dtype=float64, chunksize=(1, 1, 1830, 1830), chunktype=numpy.ndarray>\n",
       "Coordinates:\n",
       "  * time     (time) object 2017-06-28 2017-07-01 ... 2021-02-25 2021-02-27\n",
       "  * band     (band) <U3 'B01'\n",
       "  * y        (y) float64 4.5e+06 4.5e+06 4.5e+06 ... 4.39e+06 4.39e+06 4.39e+06\n",
       "  * x        (x) float64 5e+05 5.001e+05 5.001e+05 ... 6.097e+05 6.098e+05"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xr.set_options(display_style='text')\n",
    "da = stack_stac(stac_items,'B01')\n",
    "da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "extended-munich",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:py_geo]",
   "language": "python",
   "name": "conda-env-py_geo-py"
  },
  "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}