Last active
December 14, 2020 07:35
-
-
Save rmg55/1acf804ef1af0c7934b265b3a653a486 to your computer and use it in GitHub Desktop.
intake_stac + stac projection extension + gdal vrt + gdal buildvrt to combine intake stac item collection to a single object
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stderr", | |
| "output_type": "stream", | |
| "text": [ | |
| "/opt/conda/envs/py_geo/lib/python3.7/site-packages/intake/source/discovery.py:138: FutureWarning: The drivers ['sql_cat', 'sql', 'sql_auto', 'sql_manual'] do not specify entry_points and were only discovered via a package scan. This may break in a future release of intake. The packages should be updated.\n", | |
| " FutureWarning)\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "import dask\n", | |
| "dask.config.set({'distributed.dashboard.link':'https://localhost:8888/proxy/{port}/status'})\n", | |
| "import intake\n", | |
| "import jinja2 as jj2\n", | |
| "import satsearch\n", | |
| "from rasterio.transform import Affine\n", | |
| "from rasterio.crs import CRS\n", | |
| "import xarray as xr\n", | |
| "xr.set_options(display_style='text')\n", | |
| "from osgeo import gdal\n", | |
| "from tempfile import NamedTemporaryFile\n", | |
| "from dask.distributed import LocalCluster,Client\n", | |
| "import time" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/html": [ | |
| "<table style=\"border: 2px solid white;\">\n", | |
| "<tr>\n", | |
| "<td style=\"vertical-align: top; border: 0px solid white\">\n", | |
| "<h3 style=\"text-align: left;\">Client</h3>\n", | |
| "<ul style=\"text-align: left; list-style: none; margin: 0; padding: 0;\">\n", | |
| " <li><b>Scheduler: </b>tcp://127.0.0.1:34965</li>\n", | |
| " <li><b>Dashboard: </b><a href='https://localhost:8888/proxy/8787/status' target='_blank'>https://localhost:8888/proxy/8787/status</a></li>\n", | |
| "</ul>\n", | |
| "</td>\n", | |
| "<td style=\"vertical-align: top; border: 0px solid white\">\n", | |
| "<h3 style=\"text-align: left;\">Cluster</h3>\n", | |
| "<ul style=\"text-align: left; list-style:none; margin: 0; padding: 0;\">\n", | |
| " <li><b>Workers: </b>12</li>\n", | |
| " <li><b>Cores: </b>12</li>\n", | |
| " <li><b>Memory: </b>26.70 GB</li>\n", | |
| "</ul>\n", | |
| "</td>\n", | |
| "</tr>\n", | |
| "</table>" | |
| ], | |
| "text/plain": [ | |
| "<Client: 'tcp://127.0.0.1:34965' processes=12 threads=12, memory=26.70 GB>" | |
| ] | |
| }, | |
| "execution_count": 2, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "clust = LocalCluster(threads_per_worker=1)\n", | |
| "cl = Client(clust)\n", | |
| "cl" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "## Stac Item Collection from SatSearch\n", | |
| "\n", | |
| "URL='https://earth-search.aws.element84.com/v0'\n", | |
| "results = satsearch.Search.search(url=URL,\n", | |
| " collections=['sentinel-s2-l2a-cogs'],\n", | |
| " bbox=[-104.79107047, 40.78311181, -104.67687336, 40.87008987],\n", | |
| " query={'eo:cloud_cover': {'lt': 10}})\n", | |
| "catalog = intake.open_stac_item_collection(results.items())" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 4, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "Starting\n", | |
| "STAC collection to VRT files completed, time:5.89 seconds\n", | |
| "Stacking VRT Files in to single Dataset, time:0.1 seconds\n", | |
| "Lazy Load of data into DataArray, time:0.02 seconds\n", | |
| "\n" | |
| ] | |
| }, | |
| { | |
| "data": { | |
| "text/html": [ | |
| "<pre><xarray.DataArray (obs_date: 169, y: 10980, x: 10980)>\n", | |
| "dask.array<open_rasterio-e9a1f8c7d57c3e928db71965155dedb6<this-array>, shape=(169, 10980, 10980), dtype=uint16, chunksize=(1, 5490, 5490), chunktype=numpy.ndarray>\n", | |
| "Coordinates:\n", | |
| " * obs_date (obs_date) object 2020-10-13T18:02:42+00:00 ... 2017-01-19T17:4...\n", | |
| " * y (y) float64 4.6e+06 4.6e+06 4.6e+06 ... 4.49e+06 4.49e+06 4.49e+06\n", | |
| " * x (x) float64 5e+05 5e+05 5e+05 ... 6.098e+05 6.098e+05 6.098e+05\n", | |
| "Attributes:\n", | |
| " transform: (10.0, 0.0, 499980.0, 0.0, -10.0, 4600020.0)\n", | |
| " crs: +init=epsg:32613\n", | |
| " res: (10.0, 10.0)\n", | |
| " is_tiled: 1\n", | |
| " nodatavals: (nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,...\n", | |
| " scales: (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,...\n", | |
| " offsets: (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,...</pre>" | |
| ], | |
| "text/plain": [ | |
| "<xarray.DataArray (obs_date: 169, y: 10980, x: 10980)>\n", | |
| "dask.array<open_rasterio-e9a1f8c7d57c3e928db71965155dedb6<this-array>, shape=(169, 10980, 10980), dtype=uint16, chunksize=(1, 5490, 5490), chunktype=numpy.ndarray>\n", | |
| "Coordinates:\n", | |
| " * obs_date (obs_date) object 2020-10-13T18:02:42+00:00 ... 2017-01-19T17:4...\n", | |
| " * y (y) float64 4.6e+06 4.6e+06 4.6e+06 ... 4.49e+06 4.49e+06 4.49e+06\n", | |
| " * x (x) float64 5e+05 5e+05 5e+05 ... 6.098e+05 6.098e+05 6.098e+05\n", | |
| "Attributes:\n", | |
| " transform: (10.0, 0.0, 499980.0, 0.0, -10.0, 4600020.0)\n", | |
| " crs: +init=epsg:32613\n", | |
| " res: (10.0, 10.0)\n", | |
| " is_tiled: 1\n", | |
| " nodatavals: (nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,...\n", | |
| " scales: (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,...\n", | |
| " offsets: (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,..." | |
| ] | |
| }, | |
| "execution_count": 4, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "## Stack Bands across items by:\n", | |
| "## 1. Converting StacItems to VRT objects\n", | |
| "## 2. band-stack VRTS\n", | |
| "## 3. load into xarray DataArray.\n", | |
| "\n", | |
| "\n", | |
| "##### Build VRT Template #####\n", | |
| "vrt_template = jj2.Template('''\n", | |
| "<VRTDataset rasterXSize=\"{{rasterXSize}}\" rasterYSize=\"{{rasterYSize}}\">\n", | |
| " <SRS>{{SRS}}</SRS>\n", | |
| " <GeoTransform>{{GeoTransform}}</GeoTransform>\n", | |
| " <VRTRasterBand dataType=\"{{dtype}}\" band=\"1\">\n", | |
| " <SimpleSource>\n", | |
| " <SourceFilename relativeToVRT=\"1\">/vsicurl/{{SourceFilename}}</SourceFilename>\n", | |
| " <SourceBand>1</SourceBand>\n", | |
| " <SourceProperties RasterXSize=\"{{rasterXSize}}\" RasterYSize=\"{{rasterYSize}}\" DataType=\"{{dtype}}\" BlockXSize=\"1024\" BlockYSize=\"1024\" />\n", | |
| " <SrcRect xOff=\"0\" yOff=\"0\" xSize=\"{{rasterXSize}}\" ySize=\"{{rasterYSize}}\" />\n", | |
| " <DstRect xOff=\"0\" yOff=\"0\" xSize=\"{{rasterXSize}}\" ySize=\"{{rasterYSize}}\" />\n", | |
| " </SimpleSource>\n", | |
| " </VRTRasterBand>\n", | |
| "</VRTDataset>\n", | |
| "''')\n", | |
| "\n", | |
| "#### Loop over Items in catalog and build a VRT representation ####\n", | |
| "Band='B04'\n", | |
| "vrt_lst = []\n", | |
| "band_mapping = {}\n", | |
| "print(\"Starting\")\n", | |
| "t1 = time.time()\n", | |
| "for i,item_name in enumerate(list(catalog)):\n", | |
| " item = catalog[item_name]\n", | |
| " if i == 0:\n", | |
| " dtype = str(item[Band].to_dask().dtype)\n", | |
| " band_mapping[i]=item.metadata['datetime']\n", | |
| " vrt = vrt_template.render(rasterXSize = item[Band].metadata['proj:shape'][0],\n", | |
| " rasterYSize = item[Band].metadata['proj:shape'][1],\n", | |
| " SourceFilename = item[Band].metadata['href'],\n", | |
| " SRS=CRS.from_epsg(item.metadata['proj:epsg']).wkt,\n", | |
| " GeoTransform=', '.join(str(x) for x in Affine.to_gdal(Affine(*item[Band].metadata['proj:transform'][0:6]))),\n", | |
| " band=Band,\n", | |
| " obs_date=item.metadata['datetime'],dtype=dtype)\n", | |
| " vrt_lst.append(vrt)\n", | |
| "t2 = time.time()\n", | |
| "print('STAC collection to VRT files completed, time:'+str(round(t2-t1,2))+' seconds')\n", | |
| "\n", | |
| "\n", | |
| "#### Stack the vrt objects via gdalbuildvrt ####\n", | |
| "with NamedTemporaryFile() as tmpfile:\n", | |
| " vrt_options = gdal.BuildVRTOptions(separate=True,bandList=[1])\n", | |
| " my_vrt = gdal.BuildVRT(tmpfile.name, vrt_lst, options=vrt_options)\n", | |
| " my_vrt = None\n", | |
| " f = tmpfile.read().decode(\"utf-8\")\n", | |
| "t3 = time.time()\n", | |
| "print('Stacking VRT Files in to single Dataset, time:'+str(round(t3-t2,2))+' seconds')\n", | |
| "\n", | |
| "#### Open vrt file in Xarray ####\n", | |
| "ds = xr.open_rasterio(f,chunks={'band':1,'x':'auto','y':'auto'})\n", | |
| "ds = ds.rename({'band':'obs_date'})\n", | |
| "ds['obs_date'] = [band_mapping[b-1] for b in ds.obs_date.values.tolist()]\n", | |
| "print('Lazy Load of data into DataArray, time:'+str(round(time.time()-t3,2))+' seconds\\n')\n", | |
| "ds" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 5, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "[<matplotlib.lines.Line2D at 0x7fefddb12090>]" | |
| ] | |
| }, | |
| "execution_count": 5, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEHCAYAAABCwJb2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAd8UlEQVR4nO3df5wc9X3f8dfn7iQhjAkgCSpL2KINDQWnwdFF5WKXx8WygTZ5GDUOrnjYQbFJVGHsoLiPUGibPtzysERcP1zF8cP4FIyRGmqsgPMAxyEYX3yPONYFegICBkzBgQoZDDI1NXZj/Tg+/eP73cfOzc3uzt7O7s7evJ+Pxzx2d3Zm9juzs5/5zuf73Rlzd0REpBqG+l0AERHpHQV9EZEKUdAXEakQBX0RkQpR0BcRqZCRfheglZUrV/q6dev6XQwRkYGxcuVK7r333nvd/ZL0e6UP+uvWrWNmZqbfxRARGShmtjJrvNI7IiIVoqAvIlIhCvoiIhWioC8iUiEK+iIiFaKgLyJSIQr6IgNuehp27gyPIq2Uvp++iDQ2PQ0bN8LRo7B0KUxOwthYv0slZaaavsgAm5oKAX92NjxOTfW7RFJ2CvoiA2x8PNTwh4fD4/h4v0skZaf0jsgAGxsLKZ2pqRDwldqRVhT0RQbc2JiCveSn9I6ISIUo6IuIVIiCvohIhSjoi4hUiIK+iEiFKOiLiFSIgr6ISIUo6IuIVIiCvohIhSjoi4hUiIK+iEiFKOiLiFSIgr6ISIUo6IuIVIiCvohIheQK+mb2O2b2mJl9y8y+YGYnmNlpZnafmT0VH09NTH+9mT1tZk+a2cWJ8evN7NH43qfMzLqxUiIikq1l0DezNcBvA6Pu/mZgGNgMXAdMuvvZwGR8jZmdG98/D7gE+IyZDcfF3QRsBc6OwyWFro2IiDSVN70zAiw3sxHgROB54FJgT3x/D7ApPr8UuN3dj7j7M8DTwAYzWw2c7O7T7u7A3sQ8IiLSAy2Dvrt/F/gEcBB4Afi/7v5V4Ax3fyFO8wJwepxlDfBcYhGH4rg18Xl6/DxmttXMZsxs5vDhw+2tkYiINJQnvXMqofZ+FvAG4HVm9r5ms2SM8ybj54903+3uo+4+umrVqlZFFBGRnPKkd94BPOPuh939GPAl4BeBF2PKhvj4Upz+EHBmYv61hHTQofg8PV5ERHokT9A/CFxgZifG3jYbgSeAu4EtcZotwF3x+d3AZjNbZmZnERpsH4gpoFfN7IK4nCsS84iISA+MtJrA3e83szuAB4HjwEPAbuAkYJ+ZXUk4MFwWp3/MzPYBj8fpr3b32bi4q4BbgeXAPXEQEZEesdCRprxGR0d9Zmam38UQERkoZnbA3UfT4/WPXBGRClHQFxGpEAV9EZEKUdAXkb6ZnoadO8Oj9EbL3jsiIt0wPQ0bN8LRo7B0KUxOwthY/b2pKRgfr4+TYijoi0hfTE2FgD87Gx6npkKAb3YwkM4pvSMifTE+HoL68HB4HB8P47MOBlIc1fRFpC/GxkItPp3GqR0MajX92sFAiqGgLyJ9MzY2P3XT6GAgxVDQF5HSyToYSDGU0xcR6bF+dlVVTV9EpIf63TtJNX0RkR7qd+8kBX0RkR5q1FW1V5TeERHpoX73TlLQl7bpL/Iineln7yQFfWlLvxuhRKQzyulLW/rdCCUinVHQl7b0uxFKRDqj9I60pd+NUCLSGQV9aZv+Ii8yuJTeERGpEAV9EZEKUdAXEakQBX0RkQpR0BcRqRAFfRGRClHQFxGpEAV9EZEKUdAXEakQBX0RkQpR0BeRrurnTcBlPl17R0S6RvdfKB/V9EWka3T/hfJR0BeRrtH9F8pH6R0R6Rrdf6F8FPRFpKt0/4VyUXpHRKRCFPRFRCokV9A3s1PM7A4z+7aZPWFmY2Z2mpndZ2ZPxcdTE9Nfb2ZPm9mTZnZxYvx6M3s0vvcpM7NurJSIiGTLW9P/A+Av3P0c4OeAJ4DrgEl3PxuYjK8xs3OBzcB5wCXAZ8xsOC7nJmArcHYcLiloPUREJIeWQd/MTgYuBD4H4O5H3f0V4FJgT5xsD7ApPr8UuN3dj7j7M8DTwAYzWw2c7O7T7u7A3sQ8IiLSA3lq+v8QOAx83sweMrObzex1wBnu/gJAfDw9Tr8GeC4x/6E4bk18nh4vIiI9kifojwA/D9zk7m8BfkxM5TSQlaf3JuPnL8Bsq5nNmNnM4cOHcxRRRETyyBP0DwGH3P3++PoOwkHgxZiyIT6+lJj+zMT8a4Hn4/i1GePncffd7j7q7qOrVq3Kuy4iItJCy6Dv7t8DnjOzn4mjNgKPA3cDW+K4LcBd8fndwGYzW2ZmZxEabB+IKaBXzeyC2GvnisQ8IiLSA3n/kfth4DYzWwr8HfB+wgFjn5ldCRwELgNw98fMbB/hwHAcuNrdZ+NyrgJuBZYD98RBRER6xEJHmvIaHR31mZmZfhdDRGSgmNkBdx9Nj9c/ckVEKkRBX0SkQhT0RUQqREFfRKRCFPRFRCpEQV9EpEIU9EVEKkRBX0SkQhT0RUQqREFfRKRCFPRFRCpEQV9EpEIU9EVEKkRBX0SkQhT0RUQqREFfRKRCFPRFRCpEQV9EpEIU9EVEKkRBX0SkQhT0RUQqREFfRKRCFPRFRCpEQV9EAJiehp07w6MsXiP9LoCI9N/0NGzcCEePwtKlMDkJY2P9LpV0g2r6IsLUVAj4s7PhcWqq3yWSblHQFxHGx0MNf3g4PI6P97tE0i1K74gIY2MhpTM1FQK+UjuLl4K+iAAh0CvYL35K74iIVIiCvohIhSjoi4hUiIK+iEiFKOiLiFSIgr6ISIUo6IuIVIiCvohIhSjoi4hUiIK+iEiF5A76ZjZsZg+Z2Z/F16eZ2X1m9lR8PDUx7fVm9rSZPWlmFyfGrzezR+N7nzIzK3Z1RESkmXZq+tcATyReXwdMuvvZwGR8jZmdC2wGzgMuAT5jZsNxnpuArcDZcbiko9KLiEhbcgV9M1sL/DJwc2L0pcCe+HwPsCkx/nZ3P+LuzwBPAxvMbDVwsrtPu7sDexPziIhID+St6e8CrgVeS4w7w91fAIiPp8fxa4DnEtMdiuPWxOfp8SIi0iMtg76Z/QrwkrsfyLnMrDy9Nxmf9ZlbzWzGzGYOHz6c82NFRKSVPDX9twLvMrNngduBt5vZHwMvxpQN8fGlOP0h4MzE/GuB5+P4tRnj53H33e4+6u6jq1atamN1RESkmZZB392vd/e17r6O0ED7l+7+PuBuYEucbAtwV3x+N7DZzJaZ2VmEBtsHYgroVTO7IPbauSIxj4iI9EAnd866EdhnZlcCB4HLANz9MTPbBzwOHAeudvfZOM9VwK3AcuCeOIiISI9Y6EhTXqOjoz4zM9PvYoiIDBQzO+Duo+nx+keuiEiFKOiX1PQ07NwZHkVEitJJTl+6ZHoaNm6Eo0dh6VKYnISxsX6XSkQWA9X0S2hqKgT82dnwODXV7xKJyGKhoF9C4+Ohhj88HB7Hx/tdIhFZLJTeKaGxsZDSmZoKAV+pHREpioJ+SY2NKdiLSPGU3hERqRAFfRGRClHQFxGpEAV9EZEKUdAXEakQBX0RkQpR0BcRqRAFfRGRClHQFxGpEAV9EZEKUdAXEakQBX0RkQpR0BdJ0B3LZLHTVTZFIt2xTKpANX2RqKg7lulsYbBU7ftSTV8kqt2xrFbTX8gdy3S2MFiq+H2ppi8S1e5YdsMNC//x6/7Gg6WK35dq+iIJyTuWTU+3f8vKIs4WpHeq+H0p6ItkWOhpv+5vPFiq+H0p6ItkyDrtzxsQdH/j8so6e6va96WgL5Khiqf9i10VG22zqCFXJEMRjbpSLr1stC1zN1DV9EUaqNpp/2LXq7O3sp9RqKYv0qYy1+KKsFjXr1dnb2XvBqqavkgbyl6L69RiX79enL2VvT1INX2RNpS9Ftepxb5+vVD29iDV9EXaUPZaXKcW+/r1SpnbgxT0Rdqw2P/Ms9jXT8Dcvd9laGp0dNRnZmb6XQwRkYYWcsmObjOzA+4+mh6vmr6I9F0Zg2Zeg9b4raAvIn01aEEzrZNLdvSDeu+ISF8Neo+hWuP38PBgNH6rpi8dGeTTcimHQe8xNGiN3y2DvpmdCewF/gHwGrDb3f/AzE4DvgisA54F3uPuP4jzXA9cCcwCv+3u98bx64FbgeXAnwPXeNlbkqWhQT8tl3IYtKCZpZMumr2uOOWp6R8H/q27P2hmrwcOmNl9wG8Ak+5+o5ldB1wH/DszOxfYDJwHvAH4mpn9Y3efBW4CtgJ/Qwj6lwD3FL1S0huDlsuU8ipzv/Zu6kfFqWVO391fcPcH4/NXgSeANcClwJ442R5gU3x+KXC7ux9x92eAp4ENZrYaONndp2Ptfm9iHhlAg5bLlP5arNf0SWtnPfvRntFWTt/M1gFvAe4HznD3FyAcGMzs9DjZGkJNvuZQHHcsPk+Pz/qcrYQzAt74xje2U0TpocVwWi7dkU5ZVCUV2Gg9G6Vw+tGekTvom9lJwJ3Adnf/oZk1nDRjnDcZP3+k+25gN4Q/Z+Uto/ReVU/LpbGswLfYUoGNgnijmntye+zaBS+/XJ931y648044//y5Nf1uVaZyBX0zW0II+Le5+5fi6BfNbHWs5a8GXorjDwFnJmZfCzwfx6/NGC8ii0hW4Bv0HjpJzc5astYzuT3+/u9h2zYwg2XLQsDfvh2OHIGvfhWGhkK69LXXwD1MU/RZUcucvoUq/eeAJ9z9k4m37ga2xOdbgLsS4zeb2TIzOws4G3ggpoJeNbML4jKvSMwjsuhUJYedltXWU/YrT7ajWR4+az3Hx2EkUb12D0H9yJFQwz96NLyG8HjsWFh2bZqi8/x5avpvBX4deNTMHo7j/j1wI7DPzK4EDgKXhRXyx8xsH/A4oefP1bHnDsBV1Lts3oN67lRG1frzVyWHnaVRW89iSQXWDmpHjoSa+YoVc99Pr+fYGLz//fDZz86dbmgI3v1u+MY3wrJqgT/JrPizopZB393/mux8PMDGBvN8DPhYxvgZ4M3tFFAG32IOgO3kdhfLOudRVIAvU2UhWZZdu+Dqq8P3u307/OzPNi5f7UxvyZJQi4cQzD/yEdi6Ncy7fTs88MD8eRs3nS6c/pErXbdYA2C7uV1pT5kqC+mybNlST9MkUzzpA1RtviNHQgA//3x49NEw7x/+IWzaFKb727+tf9bQUL3W717870XX3pGu61Z//lY5827n1NvN7Q6y3bvh4ovDY69MTYVgWWsA3b69f+0je/fCT35S/66hvk+PjIRa+vg4/N7vhSA/PR2G7dtD2V97Lcz7yCNzDxZ798JHPzr3DOBd74Lly7v4/xd3L/Wwfv16l8G3f7/7jh3hsajlLV/uPjwcHtPLbfV+L8ow6Grf2Xvf6x5CVRgmJrKnK3r9Jybmfi64L13a++28f7/7smXzy7B/v/umTeH7T5ZxeNh927a589QGM/eRkTDNsmVhWUND4b2hofp+VMQ2BWY8I6YqvSMdyZtzLbIRb3o61I5qjV9ZKaNepJT69ee0XuS5k2mJdAPjnXeGXDSEmv/VV4dpiu5e+PLL88cdO9b4u+zWdpmamlsT/8AH6n+4+spXwj5WYxZq51A/I0i+d8IJ9X76Bw/CH/1R2HZDQ/COd4T9Otnw3Q0K+rJgvci5pn/ItSAzOxvqTmbZPSh6lVPvdY+U9DZP/9GnKLXUSlaPkvPPr5flQx+C48fD61r3wmQ+u5MgPD4+t/ETQsoj67tcyL6Yt3wrVszNsZ98ckgbHjw4N+APD8Nv/RZccUV4/fnPh20CYT2uvDK8l9w+e/bUy5wM+N2koC8Llq5N792b70eU98eW/iF/+MPwiU/Uf4BmYUj3oKgtPxkQIfxQy9ALJI88vYKOHAkHQPfOD7q1z1uxImyzV17JDvhQb4Ccmpob9IaG6tu6iArB2FgIlBMT9QP8b/5m9nLa3RfbKd9DD4XPrpXhk58Mz0dGwnD8eAj4n/50/QwI4OtfD+WAucE+uX59uYxJVs6nTMNizOl3Kwfaa8mcdi0/2Sq/PTERcprJ/GUjO3bU86VDQ/XcZ21Ijhserm/TdJ49PW5iotzbv1lbQfK9JUvmr/9C1L4Ts/p2XbJk7uuf/un6Z5mFnHVyviVL3K+9tr5dd+yYm6teaNnytpu0uy8m962sbbd/f1jHTZvCuiXz9cltvm1befclGuT0+x7UWw0LDfplC6y18kxMlK/xr5NtVftxbNjQOgDt3z/3B2TmftFF+X7IySBUm/fCC+vjly2rr0f6x5w+eNQOOsPD4Ufdq+8g73bOWofkvEXuS+nvJHlArTU41g6UyYbJ4eEwJAN+sizXXjt3eckDQm2f2bYtX5nzbrfadNu2NQ/otWmbHVgbNcJu2lS+328jlQr6eXp29PKA0I3aWSdlSf7g2umBkrXdavNn9UBIS9b+ksGlVQ0uGeBqtfuhobkHgVqPiqwziWYHDwjTd3tfaHc7p89MmgWorO8kfYBo9HnNvpP0GdG2bfO3XS0YJs8EhofDwTy5TyR7rCQPMt3ojdPO2UHWttmxY/56mhXbs6YXKhX0m526FdnNbiE1t3QNqhc7UTJwprue5akV1ZaRtd3S65an5p4O2nkPfrUDVrqLXPJUu7b8dC0+uQ2yarbbtrW/XbPK1+h7bJVOaLasduZtN82R/E7SKZpG0zYK/Mn9KnmgSlZ0soJpNyo+nZ69Jn8nS5bkPyspk0ZBf1E25DbruVFUV75mDUHpRrh0eXbtCo1DEP6dt317vuUsRLKcQ0Nze0IcOwYPPhgaoaB5L5dG2y29bs16ICQbrlasmLveeXrXjI2Fed2z3//e9+ZevOquu+Dee+du06kpeM974LbbWn9eO9JdF9O9atrtTZTuFdRq3tq+cvBg/XtK9jhptK+305hYuwzwBz9Yb8CtXSYg+Z2YhUb1rO96ZKR+UTEIvVq60bOqk15VY2OtG2EHWtaRoExD0Tn9omr6jWpfjZY/MRFqwRMT+dI9Cyln1jqna+LpWvLQUKjVtKrJtMqBLqRW1Wy+rLxvbdzSpdk1xnTqINmImE5DJYdae8BCy5zVVpE+m2u0TkVtq6zafd7G9XYk96dkjjv5fTRqGF1oTl8Whiqld1opIp2SJ93RqEdJMkWRle5ptJyFlCfr9D1vw2s3tlue5WeloZI5/eHh0IibDv5DQ+5r184N6LV8fXJ7Jqdvlo5Kl23btvkNx+7z8+JZvYpqy+hGI2B6X0n2KElWNpqtW97vNWsdkgfjQWjgrAoF/S7I+rGkGxJrP7p0N69WXQjbbYxqlJvfv7/+V/FGjZtl+JGmz37StfgNG+YG+JGRsN1qwSadM671KqkFu7wNzo1q8lm57Fo7wP799QPQ8PD8niwLPZAvZNtlHfCbnZ0tpAdQs7PoQWjgrAoF/R7IqlU3CjR5fiCtpmnVaJcVrNIHhLL8SJM9Q9JpqCVLGqdt3EP5071FsmrwyUDXTuov6ywhHfSXLZt7BtDs4NGNA22r1F6jM47kevWjN5l0j4J+D6R/ZBddlL9nS1LeYNzstD79frrbWZmke0vU0jm1vG+6x07WehQRUFsFyWQKJ9nVsN3eNb060LZ7EOvHxcykexT0eyD9I1voqfNC+3M3q7nmaaztl3TjYLr7ZN716DSg5k2H9LIG36m86aqs7S6DrVHQt/BeeY2OjvrMzEy/i5Fbuptlu90ud+4M1+SenQ1dKW+4Aa6/Pv/ntft+GeS5iFiv1mOhnzMI2zlpejp0SbzllrCv9fsmJVI8Mzvg7qPzxivol0uZ7hbUS8kLfjX634IUb9AOVpJfo6C/KP+cNcj6duW9Pqv9mWbnzsV5a8WyWiw3K5f8FPRLqMo/RN1bVqS7FPSlVKp6piPSKwr6UjpVPtMR6bahfhdApB+mp0P7wfR0v0si0luq6UvlVLWHlAiopi8VlHWZaJGqUNCXyqn1EBoeVg8hqR6ld6Ry1ENIqkxBXypJPYSkqpTeERGpEAV9EZEKUdAXEakQBX0RkQpR0BcRqRAFfRGRCin9TVTM7DDwv3NOvhL4fheL06mylw9UxqKUvYxlLx+ojJ34PoC7X5J+o/RBvx1mNpN1p5iyKHv5QGUsStnLWPbygcrYLUrviIhUiIK+iEiFLLagv7vfBWih7OUDlbEoZS9j2csHKmNXLKqcvoiINLfYavoiItKEgr6ISJW4e98G4Ezg68ATwGPANXH8acB9wFPx8dQ4fkWc/kfApxPLeT3wcGL4PrCrwWeuBx4FngY+RT3FdSHwIHAc+LUSlu+/Jeb/X8ArRZYxvnd5/OxHgL8AVpZlGxZUvsxt2IVy/utYxseAjzfZ//u1HTstX1H74juBA/EzDgBvb/XZnWzDEpax4f7Y1bjbiw9pslOtBn4+Pn99XPFzgY8D18Xx1wG/H5+/DngbsC29I6eWewC4sMF7DwBjgAH3AP8ijl8H/FNgL/UfWmnKl5rmw8AtRZaRcG+Fl4iBNM7/0bJswyLK12gbFlzOFcBBYFV8vQfYWKLt2HH5CtwX3wK8IT5/M/Dddj673W1YtjI22x+7OfQ16Ges+F2EI+uTwOrEl/RkarrfoEFQBc4GniPjqBuX9e3E68uBidQ0tyZ3krKVL47fD7yzyDICS4DDwJvizvlZYGtZtmGR5Wu1DTss5y8AX0u8/nXgMyXajoWVr6h9MY434GVgWRuf3dE2LEsZ8+yPRQ6lyemb2TrCUfV+4Ax3fwEgPp7exqIuB77ocUumrAEOJV4fiuMGpnxm9ibgLOAviyyjux8DriKchj5PqP18biFlbKQs5Wu2DTstJ+EU/hwzW2dmI8AmQkqh7XI2UpbyFbwvvht4yN2P5PnsvGVspixlbLU/Fq0UQd/MTgLuBLa7+w87XNxm4AuNPipjXFbwnTtTucq3GbjD3WeLLKOZLSEE1bcAbyDkfK9fYBmzll+m8mVuwyLK6e4/iOX8IvAN4FlCXnkh5Zw/U7nKV8i+aGbnAb8P/Js2Prud6bI+s0xlbLg/dkPfg378Md8J3ObuX4qjXzSz1fH91YRcbp5l/Rww4u4H4uthM3s4Dv+FcJRdm5hlLaHWOEjlm3fQKKiM5wO4+3fiWcg+4BdLtA2LLF/mgbeo79rdv+zu/8zdxwhpg6dKtB2LLF/H+6KZrQX+FLjC3b8TR2d+dhHbsKRlbFYRLFxfg76ZGeEU/Ql3/2TirbuBLfH5FkLeLY/LSWw8d5919/Pj8J/iadurZnZB/Owrmi27bOUzs58BTgWmu1DG7wLnmtmq+PqdcZll2YaFlC9rGxZcTszs9Ph4KvBB4OYSbcdCylfEvmhmpwBfAa5392/WJm702Z1uwzKWsdH+2FXtNgIUORB6FzjhVP3hOPxLQg+DSUL3qUngtMQ8zwL/h9AV7RBwbuK9vwPOafGZo8C3gO8An6befeoX4vJ+TGiseaxM5YvvfRS4sVvbkNDT44m4rC8DK8q0DTstX6Nt2IXt+AXg8ThsLtu+2Gn5itoXgf8Y1/HhxHB6q89e6DYsWxmb7Y/dHHQZBhGRCul7Tl9ERHpHQV9EpEIU9EVEKkRBX0SkQhT0RUQqREFfRKRCFPSlkixce+ZbBS/zRy3eP8XMPljkZ4q0S0FfpHdOIfwDVqRvFPSlEszsI2b2rThsj6NHzGyPmT1iZneY2Ylx2hvN7PE4/hNNlnmWmU2b2f80sxsS408ys0kze9DMHjWzS+NbNwL/KF6X5b/GaX83zv+Imf3n7qy9SJ3+kSuLnpmtJ1xX/QLCVQ/vB95HuLPS29z9m2Z2C+GyBLcQroNyjru7mZ3i7q80WO7dhKsj7jWzqwk33jjJwiWLT3T3H5rZSuBvCPdReBPwZ+7+5jj/RcCvEa7caITrv3zc3f+qKxtCBNX0pRreBvypu//Y3X8EfAn458BzXr+I1h/H6X4I/AS42cx+Ffh/TZb7VuoX0PvvifEG7DCzR4CvEa6ffkbG/BfF4SHCAegcwsFBpGtG+l0AkR7IuqY5zL+uubv7cTPbAGwkXPL2Q8Dbmyw761T5vcAqYL27HzOzZ4ETGpRrp7tPNCu8SJFU05cq+Ctgk5mdaGavA/4V4QYibzSzsTjN5cBfW7i5xk+5+58D24nX8W/gm4QDA4RAX/NTwEsx4P8SIa0D8Crhvqw19wIfiJ+Jma2pXfZYpFtU05dFz90fNLNbCTeoBrgZ+AHhMs1bzGyCcEndmwgB+y4zO4FQE/+dJou+BvgfZnYN4aYcNbcBXzazGcKleL8dy/GymX0zdhW9x91/18z+CTAdLrXOjwhtDbluyiOyEGrIFRGpEKV3REQqROkdkRbM7D8Al6VG/4m7f6wf5RHphNI7IiIVovSOiEiFKOiLiFSIgr6ISIUo6IuIVMj/Byvyr1RGobOaAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "ds.sel(x=slice(517617.2187,527253.4091),y=slice(4524372.5,4514729.5)).mean(dim=['x','y']).compute().plot.line('b.')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 6, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "#### BuildVRT with Jinja rather than GDAL - Work in progress ####\n", | |
| "\n", | |
| "# template_str = \\\n", | |
| "# '''\n", | |
| "# {%for item_name in catalog%}\n", | |
| "# <VRTDataset rasterXSize=\"{{catalog[item_name][band].metadata['proj:shape'][0]}}\" rasterYSize=\"{{catalog[item_name][band].metadata['proj:shape'][1]}}\">\n", | |
| "# <SRS>{{SRS}}</SRS>\n", | |
| "# <Metadata>\n", | |
| "# <MDI key=\"obs_date\">{{catalog[item_name].metadata['date']}}</MDI>\n", | |
| "# </Metadata>\n", | |
| "# <GeoTransform>{{GeoTransform}}</GeoTransform>\n", | |
| "# <VRTRasterBand band=\"1\">\n", | |
| "# <SimpleSource>\n", | |
| "# <SourceFilename relativeToVRT=\"1\">/vsicurl/{{catalog[item_name][band].metadata['href']}}</SourceFilename>\n", | |
| "# <SourceBand>1</SourceBand>\n", | |
| "# <SourceProperties RasterXSize=\"{{rasterXSize}}\" RasterYSize=\"{{catalog[item_name][band].metadata['proj:shape'][1]}}\"/>\n", | |
| "# <SrcRect xOff=\"0\" yOff=\"0\" xSize=\"{{catalog[item_name][band].metadata['proj:shape'][0]}}\" ySize=\"{{catalog[item_name][band].metadata['proj:shape'][1]}}\" />\n", | |
| "# <DstRect xOff=\"0\" yOff=\"0\" xSize=\"{{catalog[item_name][band].metadata['proj:shape'][0]}}\" ySize=\"{{catalog[item_name][band].metadata['proj:shape'][1]}}\" />\n", | |
| "# </SimpleSource>\n", | |
| "# </VRTRasterBand>\n", | |
| "# </VRTDataset>\n", | |
| "# {% endfor %}\n", | |
| "# '''\n", | |
| "# Band='B04'\n", | |
| "# r = jj2.Template(template_str,trim_blocks=True).render(catalog=catalog,\n", | |
| "# band=Band,\n", | |
| "# SRS = CRS.from_epsg(catalog[list(catalog)[0]].metadata['proj:epsg']).wkt,\n", | |
| "# GeoTransform=', '.join(str(x) for x in Affine.to_gdal(Affine(*catalog[list(catalog)[0]][Band].metadata['proj:transform'][0:6]))),)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 7, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "\n", | |
| "<VRTDataset rasterXSize=\"10980\" rasterYSize=\"10980\">\n", | |
| " <SRS>PROJCS[\"WGS 84 / UTM zone 13N\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-105],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32613\"]]</SRS>\n", | |
| " <GeoTransform>499980.0, 10.0, 0.0, 4600020.0, 0.0, -10.0</GeoTransform>\n", | |
| " <VRTRasterBand dataType=\"uint16\" band=\"1\">\n", | |
| " <SimpleSource>\n", | |
| " <SourceFilename relativeToVRT=\"1\">/vsicurl/https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/13/T/EF/2017/1/S2A_13TEF_20170119_0_L2A/B04.tif</SourceFilename>\n", | |
| " <SourceBand>1</SourceBand>\n", | |
| " <SourceProperties RasterXSize=\"10980\" RasterYSize=\"10980\" DataType=\"uint16\" BlockXSize=\"1024\" BlockYSize=\"1024\" />\n", | |
| " <SrcRect xOff=\"0\" yOff=\"0\" xSize=\"10980\" ySize=\"10980\" />\n", | |
| " <DstRect xOff=\"0\" yOff=\"0\" xSize=\"10980\" ySize=\"10980\" />\n", | |
| " </SimpleSource>\n", | |
| " </VRTRasterBand>\n", | |
| "</VRTDataset>\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "# Show VRT representation of a STAC Item\n", | |
| "print(vrt)" | |
| ] | |
| } | |
| ], | |
| "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.7.8" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 4 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment