notebook-demos/reading-in-sentinel-data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reading in Sentinel-2 Images\n",
    "\n",
    "Sentinel-2 is an observation mission developed by the European Space Agency to monitor the surface of the Earth [official website](http://www.esa.int/Our_Activities/Observing_the_Earth/Copernicus/Sentinel-2).\n",
    "Sets of images are taken of the surface where each image corresponds to a specific wavelength. These images can provide useful data for a wide variety of industries, however, the format they are stored in can prove difficult to work with. This being `JPEG 2000` (file extension `.jp2`), an image compression\n",
    "format for JPEGs that allows for improved quality and compression ratio.\n",
    "\n",
    "## Why Use GeoPySpark\n",
    "\n",
    "There are few libraries and/or applications that can work with `jp2`s and big data, which can make processing large amounts of sentinel data difficult. However, by using GeoPySpark in conjunction with the tools available in Python, we are able to read in and work with large sets of sentinel imagery."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting the Data\n",
    "\n",
    "Before we can start this tutorial, we will need to get the sentinel images. All sentinel data can be found on Amazon's S3 service, and we will be downloading it straight from there.\n",
    "\n",
    "We will download three different `jp2`s that represent the same area and time in different wavelengths: Aerosol detection (443 nm), Water vapor (945 nm), and Cirrus (1375 nm). These bands are chosen because they are all in the same 60m resolution. The tiles we will be working with cover the eastern coast of Corsica taken on January 4th, 2017.\n",
    "\n",
    "For more information on the way the data is stored on S3, please see this [link](http://sentinel-pds.s3-website.eu-central-1.amazonaws.com/). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!curl -o /tmp/B01.jp2 http://sentinel-s2-l1c.s3.amazonaws.com/tiles/32/T/NM/2017/1/4/0/B01.jp2\n",
    "!curl -o /tmp/B09.jp2 http://sentinel-s2-l1c.s3.amazonaws.com/tiles/32/T/NM/2017/1/4/0/B09.jp2\n",
    "!curl -o /tmp/B10.jp2 http://sentinel-s2-l1c.s3.amazonaws.com/tiles/32/T/NM/2017/1/4/0/B10.jp2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Code\n",
    "\n",
    "Now that we have the files, we can begin to read them into GeoPySpark."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import rasterio\n",
    "import geopyspark as gps\n",
    "import numpy as np\n",
    "\n",
    "from pyspark import SparkContext"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "conf = gps.geopyspark_conf(master=\"local[*]\", appName=\"sentinel-ingest-example\")\n",
    "pysc = SparkContext(conf=conf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading in the JPEG 2000's\n",
    "\n",
    "`rasterio`, being backed by GDAL, allows us to read in the `jp2`s. Once they are read in, we will then combine the three seperate numpy arrays into one. This combined array represents a single, multiband raster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "jp2s = [\"/tmp/B01.jp2\", \"/tmp/B09.jp2\", \"/tmp/B10.jp2\"]\n",
    "arrs = []\n",
    "\n",
    "for jp2 in jp2s:\n",
    "    with rasterio.open(jp2) as f:\n",
    "        arrs.append(f.read(1))\n",
    "\n",
    "data = np.array(arrs, dtype=arrs[0].dtype)\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating the RDD\n",
    "\n",
    "With our raster data in hand, we can how begin the creation of a Python `RDD`. Please see the [core concepts](core-concepts.ipynb) guide for more information on what the following instances represent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create an Extent instance from rasterio's bounds\n",
    "extent = gps.Extent(*f.bounds)\n",
    "\n",
    "# The EPSG code can also be obtained from the information read in via rasterio\n",
    "projected_extent = gps.ProjectedExtent(extent=extent, epsg=int(f.crs.to_dict()['init'][5:]))\n",
    "projected_extent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may have noticed in the above code that we did something weird to get the `CRS` from the rasterio file. This had to be done because the way rasterio formats the projection of the read in rasters is not compatible with how GeoPySpark expects the `CRS` to be in. Thus, we had to do a bit of extra work to get it into the correct state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Projection information from the rasterio file\n",
    "f.crs.to_dict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The projection information formatted to work with GeoPySpark\n",
    "int(f.crs.to_dict()['init'][5:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# We can create a Tile instance from our multiband, raster array and the nodata value from rasterio\n",
    "tile = gps.Tile.from_numpy_array(numpy_array=data, no_data_value=f.nodata)\n",
    "tile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Now that we have our ProjectedExtent and Tile, we can create our RDD from them\n",
    "rdd = pysc.parallelize([(projected_extent, tile)])\n",
    "rdd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating the Layer\n",
    "\n",
    "From the `RDD`, we can now create a `RasterLayer` using the `from_numpy_rdd` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# While there is a time component to the data, this was ignored for this tutorial and instead the focus is just\n",
    "# on the spatial information. Thus, we have a LayerType of SPATIAL.\n",
    "raster_layer = gps.RasterLayer.from_numpy_rdd(layer_type=gps.LayerType.SPATIAL, numpy_rdd=rdd)\n",
    "raster_layer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Where to Go From Here\n",
    "\n",
    "By creating a `RasterLayer`, we can now work with and analyze the data within it. If you wish to know more about these operations, please see the following guides: [Layers Guide](layers.ipynb), [map-algebra-guide], [visulation-guide], and the [catalog-guide]."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.5.1+"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}