Preprocessing of Sentinel-1 SAR data via Snappy Python module

Preprocessing of Sentinel-1 SAR data via Snappy Python module

This chapter demonstrates the Snappy Python module for the automatization of the ESA SNAP tool.

Code examples will be shown for an automated processing chain for the preprocessing of Sentinel-1 SAR data including Calibration, Subsetting and Terrain Correction of GRD (Ground Range Detected data).

A detailed installation tutorial for snappy can be found here: https://senbox.atlassian.net/wiki/display/SNAP/How+to+use+the+SNAP+API+from+Python

First, import the needed Python modules:

 

Now loop through all Sentinel-1 data sub folders that are located within a super folder (of course, make sure, that the data is already unzipped):

Then, read in the Sentinel-1 data product:

If polarization bands are available, spolit up your code to process VH and VV intensity data separately. The first step is the calibration procedure by transforming the DN values to Sigma Naught respectively. You can specify the parameters to output the Image in Decibels as well.

Next, specify a subset AOI to reduce the data amount and processing time. The AOI specified by its outer polygon corners and is formatted through a Well Known Text (WKT).

Apply a Range Doppler Terrain Correction to correct for layover and foreshortening effects, by using the SRTM 3 arcsecond product (90m) that is downloaded automatically. You could also specify an own DEM product with a higher spatial resolution from a local path:

Fergana_Sentinel

Image classification with Python

Image classification with Python

We are going to classify a multitemporal image stack of MODIS NDVI time series (MOD13Q1). The stack consists of 23 bands (16-day composites) with a spatial resolution of 231m in sinusoidal projection. We want to classify the different land use types, especially to discriminate different crop types.

Install Python and required image processing and scientific programming packages:

For beginners, the basic Python IDLE is sufficient for scripting. However, another IDE, LiClipse is highly recommended: http://www.liclipse.com/

We assume that you already have created a bunch of training samples in 8bit TIFF format with distinct class labeling (1,2,3,4, etc). It is necessary that these single pixels are snapped to the pixel size of the original dataset and have the same dimensions and extent.

Now write your Training script:

And now write your classification script:

Multitemporal NDVI MODIS Stack with clipped AOI:

1

Training pixels with distinct classes are created from field sampling points and snapped to the original raster extent:

2

Classification results (left image: same extent as above):

3 4

 

 

 

Extracting the central strip from LANDSAT 7 imagery

Here is a simple Python code to extract the central strip from Landsat 7 imagery (SLC-off),  that is not affected by the SLC failure. The algorithm shrinks the striping zones through a morphological filter (erosion) and creates a new shapefile AOI that extracts the desired raster extent without striping effects. The code is based on Python for ArcGIS (arcpy) – so you require a ArcGIS license.

General steps:

  1. Loop through all Landsat 7 data folders
  2. Stack bands for each image
  3. Create a mask
  4. Erode the mask by 20 pixels
  5. Convert the mask to polygon
  6. Create a minimum bounding box
  7. Clip the original raster through the bbox

 

 

Managing raster data with PostGIS and Python

Managing raster data with PostGIS and Python

PostGIS is the spatial extension of the open source database management system PostgreSQL. It helps you to manage your data (vector and raster) within a coherent geodatabase through a variety of spatial functions. Having a spatial database, the times of data clutter and messiness are over, especially when you are dealing with big data. Initially PostGIS was created to for the handling of vector data only. However, during the recent years more and more raster handling functionalities were introduced. For a complete overview of spatial raster operators, have a look at: http://postgis.net/docs/manual-2.1/RT_reference.html

 

Download and install PostgreSQL and PostGIS

Download PostgreSQL from here: http://www.postgresql.org/download/

The installer for PostgreSQL includes the PostgreSQL server, pgAdmin III; a GUI for managing and developing your databases, and StackBuilder; a package manager that can be used to download and install additional PostgreSQL applications and drivers. From the StackBuilder GUI, select Spatial Extensions and install the respective PostGIS 2.2 Bundle for PostgreSQL.

 

Create a new spatial database

In pgAdmin, create a new database (right click: New Database):

db

and the spatial extension postgis (right click on Extensions: New Extension):

pg

This will create a first table within your DB: spatial_ref_sys that contains the coordinate systems, map projections and the spatial indeces.

 

Set up Python

Python provides awesome functionality for the automated raster processing within PostGIS. Automatization is necessary especially when you deal with a lot of data and iterative processes. Python scripting is also needed as the pgAdmin GUI does not support the access of all functions.

Download Python 2.7 from here: https://www.python.org/downloads/

Psycopg2 is a Python library that accesses the objects of the PostgreSQL server and allows the execution of PostGIS commands from Python.

Download psycopg2 from here: http://www.stickpeople.com/projects/python/win-psycopg/2.6.1/psycopg2-2.6.1.win32-py2.7-pg9.4.4-release.exe

 

PostGIS scripting with Python

Import the Python libraries:

Set up input path and a loop that goes through all TIFs in the directory:

Connect to the PostgreSQL server:

Import each raster through raster2pgsql function (coordinate system epsg code is set to 32633 UTM):

Now run any PostGIS command you like. In this example we run rescale the raster to 250m spatial resolution and reproject it from UTM 33N to WGS84 (epsg code 4326). At the end, the raster may be exported locally to *.hex data format. The export is optional, we could also convert the raster to CSV or numpy array.