SENTINEL-1 Flood mapping using Snappy

SENTINEL-1 Flood mapping using Snappy

Issued February 2019 Alex McVittie

Copyright ? 2019 SkyWatch Space Applications Inc.

Prerequisites

Before beginning, you must have the snap python libraries prepared, and the two libraries pygeoif and pyshp installed (available through pip or conda). For more information on getting snappy set up, please refer to this tutorial ( for set up. This is NOT the same package that is available through pip. Also please note that it is currently only supported on Python 2.7 and 3.4. Be sure that you are running the latest version of SNAP. This tutorial is written to use the processing framework provided in SNAP 6.0. Running an older version is not tested, and there are no guarantees that it will function properly on other versions of SNAP.

At the end of this tutorial, a full list of all SNAP processing tool function names are available for future reference, as well as a comprehensive script containing all the code snippets in one runnable section.

Introduction

Synthetic Aperture Radar (SAR) is a type of sensor that captures earth observation (EO) data at a wavelength shorter than the visible light spectrum, allowing it to penetrate through clouds, and observe during both day and night, unlike an optical sensor such as LANDSAT which cannot observe the Earth's surface at night or during cloudy conditions. This gives satellites collecting SAR data a unique advantage of continuous coverage of the earth's surface, even during inclement weather.

One of the rainiest locations on Earth, the Island of Kauai in Hawaii, USA, was recently affected by extreme weather, with 685 mm of rainfall over a 24 hour period in April 2018. By using SAR data and the tools available to us through SNAP, we can see what areas were affected by flooding through a simplified classification methodology.

Purpose of tutorial

This tutorial aims to familiarize the user with common SAR processing tools available in SNAP using snappy and python, in particular:

basic product information queries terrain correction de-speckling image calibration orbit file application subsetting of data using vector file input image reclassification using band math calculations land cover downloading exporting of product objects to files.

This tutorial aims to apply these skills with some basic flood mapping analysis.

Reading products with snappy

Reading in products with Snappy is simple - you do not even need to have the data unzipped. You'll need to download two datasets for this tutorial - a boundary of Hawaii and the Sentinel 1 product. The Sentinel 1 product can be accessed through Copernicus Scihub, the full product name is S1A_IW_GRDH_1SDV_20180415T163146_20180415T163211_021480_025003_8E79. The island boundary file can be downloaded here ()

You'll need to unzip the island_boundary.zip archive, but you do not need to unzip the S1A zip archive. For this tutorial, I put this data in a new folder named data, but you can put it anywhere you wish, just be sure to update the path information in the code.

%matplotlib inline

import numpy as np import matplotlib import matplotlib.pyplot as plt import matplotlib.colors as colors import os

import snappy from snappy import Product from snappy import ProductIO from snappy import ProductUtils from snappy import WKTReader from snappy import HashMap from snappy import GPF

# For shapefiles import shapefile import pygeoif

path_to_sentinel_data = "data/S1A_IW_GRDH_1SDV_20180415T163146_20180415T163211_021480_025003_8E79.zip"

product = ProductIO.readProduct(path_to_sentinel_data)

With our ZIP file now read into snappy as a SNAP product object, we can get some basic information printed about the product.

width = product.getSceneRasterWidth() print("Width: {} px".format(width)) height = product.getSceneRasterHeight() print("Height: {} px".format(height)) name = product.getName() print("Name: {}".format(name)) band_names = product.getBandNames() print("Band names: {}".format(", ".join(band_names)))

Width: 25220 px Height: 16774 px Name: S1A_IW_GRDH_1SDV_20180415T163146_20180415T163211_021480_025003_8E79 Band names: Amplitude_VH, Intensity_VH, Amplitude_VV, Intensity_VV

Lets create a method to show a product inline. The whole product is too big to be shown in line, but we will use it later.

def plotBand(product, band, vmin, vmax):

band = product.getBand(band) w = band.getRasterWidth() h = band.getRasterHeight() print(w, h)

band_data = np.zeros(w * h, np.float32) band.readPixels(0, 0, w, h, band_data)

band_data.shape = h, w

width = 12 height = 12 plt.figure(figsize=(width, height)) imgplot = plt.imshow(band_data, cmap=plt.cm.binary, vmin=vmin, vmax=vmax)

return imgplot

Image pre-processing

Orbit file application

Before any SAR pre-processing steps occur, the product subset should be properly orthorectified to improve accuracy. To properly orthorectify the image, the orbit file is applied using the Apply-Orbit-File GPF module. SNAP is able to generate and apply a high accuracy satellite orbit file to a product.

parameters = HashMap()

GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()

parameters.put('orbitType', 'Sentinel Precise (Auto Download)') parameters.put('polyDegree', '3') parameters.put('continueOnFail', 'false')

apply_orbit_file = GPF.createProduct('Apply-Orbit-File', parameters, product)

Clipping/subsetting images

The raw image is much larger than what we need to process, so to make the processing require less resources, it is important that we first extract a subset from the image. This step could be skipped, but processing the entire scene is unnecessary and time consuming.

To clip our image, the easiest way to do so is to convert a shape boundary file into a WKT (Well Known Text) string. This is easily obtained by using a few python modules. Our shapefile of the island looks like this:

The shapefile is provided with the tutorial, and was generated from the State of Hawaii's planning department () with a few modifications (selecting out the one island we are interested in, dissolving the individual electoral boundaries into one shape file, and reprojecting to EPSG:4326). r = shapefile.Reader("data/island_boundary2.shp") g=[] for s in r.shapes():

g.append(pygeoif.geometry.as_shape(s)) m = pygeoif.MultiPoint(g) wkt = str(m.wkt).replace("MULTIPOINT", "POLYGON(") + ")" With our WKT string generated, it is time to use it to subset our data.

 ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download