-
-
Save fwrite/7ba182242ba9e3dd8c875dce5b27da48 to your computer and use it in GitHub Desktop.
GeoTiles - GEOTiff download.ipynb
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
| { | |
| "nbformat": 4, | |
| "nbformat_minor": 0, | |
| "metadata": { | |
| "colab": { | |
| "name": "GeoTiles - GEOTiff download.ipynb", | |
| "provenance": [], | |
| "authorship_tag": "ABX9TyN+Mu/QArTjrffYe7vMMjoW", | |
| "include_colab_link": true | |
| }, | |
| "kernelspec": { | |
| "name": "python3", | |
| "display_name": "Python 3" | |
| }, | |
| "language_info": { | |
| "name": "python" | |
| } | |
| }, | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "view-in-github", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "<a href=\"https://colab.research.google.com/gist/fwrite/7ba182242ba9e3dd8c875dce5b27da48/geotiles-geotiff-download.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "2RY5UlSyDTnf" | |
| }, | |
| "source": [ | |
| "# AESB2440 - Assignment 2\n", | |
| "\n", | |
| "With this tool you will download a selection of Sentinel-2 data over your 'fieldwork' area.\n", | |
| "The [Sentinel-2 constallation](https://sentinel.esa.int/web/sentinel/missions/sentinel-2) captures most of the globe every 5 days, and higher latitudes (closer to the poles) every 2-3 days.\n", | |
| "Data is extracted from the [GeoTiles.nl Sentinel-2 data cube](https://www.ncgeo.nl/index.php/nl/ncg/ncg-talent-program/ncg-talent-program-2021/item/2833-add-context-to-your-model-sentinel-2-time-series-for-everyone), an experimental service to provide easy access to the imagery recorded over the Netherlands by Sentinel-2 since 2015.\n", | |
| "\n", | |
| "In your assignment you will use data from 2021-03-31, the most recent completely cloud free day.\n", | |
| "This notebook will help you create a GeoTiff from the data cube.\n", | |
| "There is no need to unravel the workings of this script, although you are welcome to do so.\n", | |
| "In cell 4 you will be asked to enter the approximate center of your area in RD-coordinates.\n", | |
| "Update these values first, in the code or in the form, and run the cell afterwards.\n", | |
| "To run a cell, click the little play button that appears in front of each row.\n", | |
| "\n", | |
| "Feel free to contact me if you have any issues downloading the data.\n", | |
| "\n", | |
| "Adriaan van Natijne <A.L.vanNatijne@tudelft.nl>, 2021" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "3FiNeIaCMxLq" | |
| }, | |
| "source": [ | |
| "# This will install the packages necessary to execute the script.\n", | |
| "# If you are running the script locally,\n", | |
| "# you may want to do this via the Conda package manager.\n", | |
| "!pip install --upgrade zarr xarray fsspec aiohttp rasterio" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "OGCjnqzoM2cZ" | |
| }, | |
| "source": [ | |
| "import rasterio as rio, zarr, numpy as np, xarray as xr\n", | |
| "from matplotlib import pyplot as plt\n", | |
| "# When working in Google Colab, the output has to be retrieved.\n", | |
| "# Locally, the output will be stored as 'Sentinel2.tiff' in the same directory.\n", | |
| "try:\n", | |
| " from google.colab import files\n", | |
| "except ImportError:\n", | |
| " files = lambda f: print('See', f)\n", | |
| "# Progress bar\n", | |
| "from dask.diagnostics import ProgressBar\n", | |
| "ProgressBar().register()" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "7Ma_r0KYQIEq" | |
| }, | |
| "source": [ | |
| "cube = xr.open_zarr(\n", | |
| " # Request the datacube from https://geotiles.nl/cube.zarr.\n", | |
| " # The second option, 'normalize_keys', is True by default and should mannually be set to False.\n", | |
| " # 'simplecache::' is added to keep a temporary, local, copy of the data for faster repeated access.\n", | |
| " # See https://filesystem-spec.readthedocs.io/en/latest/api.html#built-in-implementations for details.\n", | |
| " store=zarr.storage.FSStore('simplecache::https://geotiles.nl/cube.zarr', normalize_keys=False),\n", | |
| " # We would like the data set of Sentinel-2, relative orbit 51.\n", | |
| " # Currently, this is the only data set available.\n", | |
| " group='/S2/R051',\n", | |
| " # Thanks to this setting, the coordinate axes of the data cube are in RD-coordinates.\n", | |
| " decode_coords=True,\n", | |
| " # Retain the original uint16 dtype of the data cube.\n", | |
| " mask_and_scale=False,\n", | |
| " # A technicallity, instead of searching for variables, a pre-generated file is used instead.\n", | |
| " consolidated=True\n", | |
| ")" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "oMiRy6NWH6WL" | |
| }, | |
| "source": [ | |
| "## Center coordinates\n", | |
| "Enter the approximate center coordinates of your area here." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "2eJ_nOBlQWH0" | |
| }, | |
| "source": [ | |
| "#@title Center coordinates\n", | |
| "#@markdown Enter the center coordinates of your area in RD coordinates.\n", | |
| "cx = 85521 #@param {type:\"integer\"}\n", | |
| "cy = 446103 #@param {type:\"integer\"}\n", | |
| "\n", | |
| "# Round coordinates to the nearest multiple of ten.\n", | |
| "cx, cy = round(cx, -1), round(cy, -1)\n", | |
| "# Create a box around it.\n", | |
| "minx, maxx = cx -5000, cx +5000\n", | |
| "miny, maxy = cy -5000, cy +5000\n", | |
| "\n", | |
| "assert minx > cube.x.min() and maxx < cube.x.max(), 'x outside Netherlands'\n", | |
| "assert miny > cube.y.min() and maxy < cube.y.max(), 'y outside Netherlands'" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "xRmjttuoRXNg" | |
| }, | |
| "source": [ | |
| "# Select the region of interest (roi) from the data.\n", | |
| "roi_cube = cube.sel(x=slice(minx, maxx-10), y=slice(maxy-10, miny), time='2021-03-31')\n", | |
| "# Bands to be exported.\n", | |
| "bands = ['B01', 'B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B8A', 'B09', 'B11', 'B12']\n", | |
| "# GeoTiff information to be exported to QGIS.\n", | |
| "profile = rio.profiles.DefaultGTiffProfile(\n", | |
| " crs=rio.crs.CRS.from_epsg(28992),\n", | |
| " dtype='uint16',\n", | |
| " count=len(bands),\n", | |
| " width=(maxx-minx)//10,\n", | |
| " height=(maxy-miny)//10,\n", | |
| " transform=rio.transform.from_origin(minx, maxy-5, 10, 10))" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "TVVzCV2AIII-" | |
| }, | |
| "source": [ | |
| "## Quick check\n", | |
| "Below you will get a preview of the area you selected.\n", | |
| "Maybe you can already recognize some features?\n", | |
| "\n", | |
| "For question 1, create a better visualization with QGIS." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "0ENP1zS_YhOC" | |
| }, | |
| "source": [ | |
| "# Stack the RGB bands in order.\n", | |
| "TCI = roi_cube[['B04', 'B03', 'B02']].to_stacked_array('rgb', ('x', 'y'))\n", | |
| "# Calculate the quantiles for scaling (2% - 98%).\n", | |
| "TCI_quantile = TCI.chunk({'x': -1, 'y': -1}).quantile((0.02, 0.98), ('x', 'y'))\n", | |
| "# Calculate a scaled version of the image.\n", | |
| "TCI_scaled = ((TCI -TCI_quantile.sel(quantile=0.02)) / \\\n", | |
| " (TCI_quantile.sel(quantile=0.98) -TCI_quantile.sel(quantile=0.02))) \\\n", | |
| " .clip(0,1).drop_vars('quantile')\n", | |
| "\n", | |
| "# Create the output map.\n", | |
| "fig, ax = plt.subplots(1,1)\n", | |
| "TCI_scaled.plot.imshow(ax=ax)\n", | |
| "ax.set_aspect('equal')" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "wOCqZHZhK7Za" | |
| }, | |
| "source": [ | |
| "## GeoTiff export\n", | |
| "\n", | |
| "Here the full image will be downloaded, converted and exported to a GeoTiff QGIS is able to work with.\n", | |
| "Your browser may offer to open the file in an image viewer.\n", | |
| "Ignore this suggestion, download the file and open it with QGIS instead." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "4vQlLgfPVDm8" | |
| }, | |
| "source": [ | |
| "with rio.open('Sentinel2.tiff', 'w', **profile) as r_out:\n", | |
| " # Write data\n", | |
| " r_out.write(roi_cube[bands].to_array())\n", | |
| " # Add band names\n", | |
| " r_out.descriptions = bands\n", | |
| " # Add copyright information\n", | |
| " r_out.update_tags(TIFFTAG_COPYRIGHT='Contains modified Copernicus Sentinel data.')\n", | |
| "\n", | |
| "files.download('Sentinel2.tiff')" | |
| ], | |
| "execution_count": null, | |
| "outputs": [] | |
| } | |
| ] | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment