Classification of Sentinel-1 Time Series

SAR data are not especially good for vegetation classification. However they have the great advantage of being independent of cloud cover. Here we investigate the use of S1 time series over a complete growth period for thematic mapping.

Since we have no ground truth data, we will use the Canada AAFC Annual Crop Inventory, which is also on the GEE archive. In particular the 2017 inventory for an area in southern Saskatchewan. This area consists of large agricultural fields, well-defined crops, and flat terrain (a big advantage for SAR measurement).

Multilook SAR image data are not normally distributed, rather they are gamma distributed. The GEE classifiers might be expected not to work so well, so we will use Tensorflow to program a more flexible neural network classifier.

First of all we grab a time series for the region of interest over the 2017 growing season (March to October):

In [2]:
%matplotlib inline
import ee
from auxil.eeWishart import omnibus
ee.Initialize()
poly = ee.Geometry({'geodesic': False, 'type': 'Polygon', 'coordinates': [[[-105.10328600000001, 50.24327999999998], 
                                                                           [-104.66649400000001, 50.24327999999998], 
                                                                           [-104.66649400000001, 50.46604134048255], 
                                                                           [-105.10328600000001, 50.46604134048255], 
                                                                           [-105.10328600000001, 50.24327999999998]]]})
poly = ee.Geometry(poly)
collection = ee.ImageCollection('COPERNICUS/S1_GRD') \
                      .filterBounds(poly) \
                      .filterDate(ee.Date('2017-03-01'), ee.Date('2017-11-01')) \
                      .filter(ee.Filter.eq('transmitterReceiverPolarisation', ['VV','VH'])) \
                      .filter(ee.Filter.eq('resolution_meters', 10)) \
                      .filter(ee.Filter.eq('instrumentMode', 'IW')) \
                      .filter(ee.Filter.eq('relativeOrbitNumber_start',5)) \
                      .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))  
collection = collection.sort('system:time_start')

collection.size().getInfo()
Out[2]:
19

Next, we transform the collection into two multiband images with linear intensities, one pre-processed with the refined Lee filter and one not:

In [6]:
from auxil.eeRL import refinedLee
import math

def get_vvvh(image):   
    ''' get 'VV' and 'VH' bands from sentinel-1 imageCollection and restore linear signal from db-values '''
    return image.select('VV','VH').multiply(ee.Image.constant(math.log(10.0)/10.0)).exp()

timeseries_rl = collection \
            .map(get_vvvh) \
            .map(refinedLee) \
            .toBands() \
            .clip(poly) \
            .float()    

timeseries_raw = collection \
            .map(get_vvvh) \
            .toBands() \
            .clip(poly) \
            .float()          

timeseries_rl.bandNames().length().getInfo()
Out[6]:
38

The class lables are conveniently obtained from the GEE archive of the Canadian AAFC Annual Crop Inventory for the year 2017, and we append them as an additional band (band 39):

In [7]:
crop2017 = ee.ImageCollection('AAFC/ACI') \
    .filter(ee.Filter.date('2017-01-01', '2017-12-01')) \
    .first() \
    .clip(timeseries_raw.geometry())\
    .float()
labeled_timeseries_rl = ee.Image.cat(timeseries_rl,crop2017)
labeled_timeseries_raw = ee.Image.cat(timeseries_raw,crop2017)

labeled_timeseries_rl.bandNames().length().getInfo()
Out[7]:
39

Now export both images to the Google drive (cloud storage would be better, but I don't have a billing account). Note that the export scale is 30m corrresponding to the AAFC/ACI resolution:

In [8]:
drexport1 = ee.batch.Export.image.toDrive(labeled_timeseries_rl,
                  description='driveExportTask', 
                  folder = 'EarthEngineImages',
                  fileNamePrefix='labeled_timeseries_rl',scale=30,maxPixels=1e10)
drexport1.start()
drexport2 = ee.batch.Export.image.toDrive(labeled_timeseries_raw,
                  description='driveExportTask', 
                  folder = 'EarthEngineImages',
                  fileNamePrefix='labeled_timeseries_raw',scale=30,maxPixels=1e10)
drexport2.start()

After downloading from the drive to a local directory, we have:

In [9]:
!ls -l imagery/regina

insgesamt 308676
-rw-rw-r-- 1 mort mort 158687709 Mär 12 15:08 labeled_timeseries_raw.tif
-rw-rw-r-- 1 mort mort 157381818 Mär 12 15:08 labeled_timeseries_rl.tif

This displays three of the (filtered) VV bands and the last (label) band:

In [48]:
%run scripts/dispms -f imagery/regina/labeled_timeseries_rl.tif  -p [21,23,25] \
-F imagery/regina/labeled_timeseries_rl.tif -E 2 -P [39,39,39]