Last active
April 19, 2024 02:12
-
-
Save mortcanty/bbdaab2835334e949a45d7c3b7bf0396 to your computer and use it in GitHub Desktop.
S1MultitempClass.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":"s1multitempclass.ipynb","provenance":[],"collapsed_sections":[],"authorship_tag":"ABX9TyOVBWDu0Uy6zca4ImD8TLDl"},"kernelspec":{"name":"python3","display_name":"Python 3"}},"cells":[{"cell_type":"markdown","metadata":{"id":"pJrMq0aNrF-u"},"source":["#Land cover classification with Sentinel-1 time series\n","## Mort Canty\n","##ZFL Bonn, March, 2021\n"]},{"cell_type":"markdown","metadata":{"id":"4xFjivyNFe7U"},"source":["## Context\n","A big advantage of SAR satellite imagery over its optical/infrared counterparts is of course weather and solar illumination independence. But for land cover classification, SAR data are relatively insensitive to vegetation/crop differences. This can be at least partly compensated by making use of time series of SAR acquisitions over a complete growth period. In this part of the course we'll examine this possibility with Sentinel-1 images collected on the Google Earthengine. There are many alternative classifiers available, such as Bayes, random forest, support vector machine, etc. We will choose to work with a neural network, essentially because of its flexibility, but also because we will continue in the next part of the course with deep learning (i.e. neural network) methods.\n","\n","### References\n","[Géron (2019) Hands-On Machine Learning ...](https://b-ok.cc/book/5341845/f49201)\n","\n","[Canty (2109) Image Anaylsis, Classification and Change Detection ...](https://www.taylorfrancis.com/books/image-analysis-classification-change-detection-remote-sensing-morton-john-canty/10.1201/9780429464348)\n","\n","[Moreira et al. (2013) A Tutorial on Synthetic Aperture Radar](https://elib.dlr.de/82313/)\n","\n","### Software\n","\n","[Colab notebooks](https://gist.github.com/mortcanty)\n","\n","[Python scripts](https://drive.google.com/drive/folders/1D2RAEu14Zcl9NzqWkSg4J7MilHGHi0BH?usp=sharing)\n","\n"]},{"cell_type":"markdown","metadata":{"id":"YtBK8zh8DsIO"},"source":["## Feed forward neural networks in a nutshell"]},{"cell_type":"markdown","metadata":{"id":"CsvlJocyhY6d"},"source":["\n","\n","\n","\n"]},{"cell_type":"markdown","metadata":{"id":"65UTsQEvSUkZ"},"source":["__Biased input vector__:\n","\n","$$\n","g^\\top = (1,g_1,\\dots,g_N)\n","$$\n","\n","__Sigmoid activation__ function (hidden layers, classical):\n","\n","$$\n","n_j(g) = {1\\over 1+ e^{-w^h_j\\cdot g}},\\quad j = 1 \\dots L.\n","$$\n","\n","__Regularized linear unit (relu) activation__ (hidden layers, more suitable for TensorFlow's autodifferentiation and for deep learning generally):\n","\n","$$\n","n_j(g) = \\max(0,w^h_j\\cdot g),\\quad j = 1 \\dots L.\n","$$\n","\n","__Softmax activation__ function (output layer):\n","\n","$$\n","m_k(n) = { e^{w^o_k\\cdot n}\\over e^{w^o_1\\cdot n}+ \\dots e^{w^o_K\\cdot n} }, \\quad k = 1\\dots K,\n","$$\n","\n","which is called the _a posteriori_ probability of class $k$ given input $g$.\n","\n","\n","This simple model is in fact a __universal classifier__. To quote from the classic text [Neural Networks for Pattern Recognition, Bishop (1995)](https://www.amazon.com/Networks-Recognition-Advanced-Econometrics-Paperback/dp/0198538642):\n","\n","... _in the context of a classification problem, networks\n","with sigmoidal non-linearities and two layers of weights can approximate any decision boundary to\n","arbitrary accuracy. ... More generally, the capability of such networks to approximate general smooth\n","functions allows them to model posterior probabilities of class membership._\n","\n","\n","__One hot label encoding__ for training example $g(\\nu)$ in class 1:\n","\n","$$\n","\\ell^{(1)}(\\nu) = (1,0,\\dots, 0)\n","$$\n","\n","In class 2:\n","\n","$$\n","\\ell^{(2)}(\\nu) = (0,1,\\dots, 0)\n","$$\n","\n","and so on.\n"," \n","__Categorical cross entropy__ (loss function) for training example $g(\\nu)$ in class $c$:\n","\n","$$\n","E(w^h,w^o) = -\\sum_k\\ell^{(c)}_k(\\nu)\\log(m_k(\\nu))\n","$$\n","\n","The loss function in minimized over batches of training examples with respect to the parameters $w^h, w^o$ using __gradient descent__ and __backpropagation__:\n","\n","$$\n","w^o \\to w^o - {\\partial E\\over \\partial w^o}\\delta,\\quad w^h \\to w^h - {\\partial E\\over \\partial w^h}\\delta\n","$$\n","\n","where $\\delta \\ll 1$ is the __learning rate__.\n","\n","When training is completed, $m^\\top=(m_1,m_2,\\dots,m_K)$ is interpreted as a __class probabilty vector__ for a given input $g$ and its class is predicted as \n","\n","$$\n","k = {\\rm argmax}(m_1,\\dots, m_K).\n","$$"]},{"cell_type":"markdown","metadata":{"id":"d9v9nkqDXhh5"},"source":["### Now that we know what a neural network classifier is, let's program one in __Python/TensorFlow__:\n","\n"," import tensorflow as tf\n"," ffn = tf.keras.models.Sequential()\n"," ffn.add(tf.keras.layers.Dense(L, activation='relu'))\n"," ffn.add(tf.keras.layers.Dense(K, activation='softmax')\n","Finished!!! \n","\n","Well, not quite, but almost. Now we'll program it properly as a Python object class:"]},{"cell_type":"code","metadata":{"id":"DF22IwbB6Vvh"},"source":["import numpy as np \n","import pandas as pd\n","import matplotlib.pyplot as plt\n","import tensorflow as tf \n"," \n","class Ffn(object): \n"," '''High-level TensorFlow (keras) FFN classifier'''\n"," def __init__(self, Ls=[10], K=10, learning_rate=0.01):\n","# setup the network architecture with the keras sequential API \n"," self._ffn = tf.keras.models.Sequential() \n","\n","# hidden layers \n"," for L in Ls:\n"," self._ffn.add(tf.keras.layers.Dense(L, activation='relu'))\n","\n","# output layer\n"," self._ffn.add(tf.keras.layers.Dense(K, activation='softmax')) \n","\n","# initialize \n"," self._history = None \n"," self._ffn.compile(\n"," optimizer=tf.keras.optimizers.SGD(learning_rate=learning_rate),\n"," metrics=['accuracy'],\n"," loss='categorical_crossentropy')\n"," \n"," def train(self, Gs, ls, epochs=10):\n"," n_split = (2*ls.shape[0])//3\n"," self._Gs_train = Gs[:n_split,:]\n"," self._Gs_valid = Gs[n_split:,:]\n"," self._ls_train = ls[:n_split,:]\n"," self._ls_valid = ls[n_split:,:]\n"," try: \n"," self._history = self._ffn.fit(self._Gs_train, self._ls_train,\n"," validation_data=(self._Gs_valid, self._ls_valid), \n"," epochs=epochs,verbose=0)\n"," return True \n"," except Exception as e:\n"," print( 'Error: %s'%e ) \n"," return None \n"," \n"," def history(self, sfn=None):\n"," pd.DataFrame(self._history.history).plot(figsize=(8,5))\n"," plt.grid(True)\n"," plt.gca().set_ylim(0,1)\n"," if sfn is not None:\n"," plt.savefig(sfn, bbox_inches='tight') \n"," plt.show() \n"," \n"," def classify(self,Gs): \n","# predict new data \n"," Ms = self._ffn.predict(Gs)\n"," cls = np.argmax(Ms, 1) + 1 \n"," return (cls, Ms)\n","\n"," def test(self, Gs, ls):\n"," m = np.shape(Gs)[0]\n"," classes, _ = self.classify(Gs)\n"," classes = np.asarray(classes, np.int16)\n"," labels = np.argmax(np.transpose(ls), axis=0) + 1\n"," misscls = np.where(classes-labels)[0]\n"," return len(misscls)/float(m) \n","\n"," def save(self, path):\n"," self._ffn.save(path)\n"," \n"," def restore(self,path):\n"," self._ffn = tf.keras.models.load_model(path)\n"," \n","# test on random data (three bands, four classes) \n","Gs = np.random.random((1000,3)) \n","ls = np.zeros((1000,4))\n","for l in ls:\n"," l[np.random.randint(0,4)]=1.0 \n","# instantiate a classidier \n","classifier = Ffn(Ls=[10,10], K=4) \n","if classifier.train(Gs, ls):\n"," classes, probabilities = classifier.classify(Gs)\n"," print( classes[:10] ) "],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"4zReZ_srIHBp"},"source":["%matplotlib inline\n","tf.keras.utils.plot_model(classifier._ffn,show_shapes=True,dpi=96)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"JbAfkmrTFajj"},"source":["## Speckle filtering\n","\n","It will turn out that the quality of the classification result is considerably improved if the images in the series are de-speckled first. So before considering the problem of classification of SAR images, we will first have a quick look at their statistical properties and discuss speckle filtering.\n","\n","#### Statistics or SAR imagery:\n","\n","https://developers.google.com/earth-engine/tutorials/community/detecting-changes-in-sentinel-1-imagery-pt-1\n","\n","### Refined Lee filter\n","\n","In the Sentinel-1 _Interferometric Wide Swath_ acquisition mode, the pixels are about 20m $\\times$ 4m (azimuth $\\times$ range) in extent and the swath widths are about 250km. For the multi-looking procedure, five cells are incoherently averaged in the range direction to achieve approx. $20m \\times 20m$ resolution.\n","\n","The observed signal $s$ for the multi-look SAR intensity image is _gamma_ distributed \n","\n","$$\n","p_{\\gamma;\\alpha,\\beta}(s) = Cs^{\\alpha-1}e^{-s/\\beta}\n","$$\n","\n","where $\\alpha, \\beta$ are parameters and $C$ is a normalization constant.\n","\n","The mean and variance are given by\n","\n","$$\n"," \\langle s\\rangle = \\mu = \\alpha\\beta,\\quad {\\rm var}(s) = \\alpha\\beta^2 = \\mu^2/\\alpha,\\quad \\alpha=m,\n","$$\n","\n","where $m$ is the _number of looks_, in our case $m=5$. \n"]},{"cell_type":"code","metadata":{"id":"TZiHPr3P03ZX"},"source":["import scipy.stats as st\n","\n","alpha = 5\n","mu = 0.1\n","beta = mu/alpha\n","z = np.linspace(0,0.5,100)\n","for alpha in range(3,7):\n"," plt.plot(z, st.gamma.pdf(z, alpha, 0, beta),label = str(alpha))\n","plt.legend() "],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"01Hat7XBXokt"},"source":["The signal can be expressed as a _multipicative speckel model_\n","$$\n","s = \\mu\\cdot v\n","$$\n","where $\\mu$ is the underlying _signal strength_ (reflectance) and $v$ represents the multiplicative speckle, with mean value $\\langle v\\rangle=1$ and variance ${\\rm var}(v)=1/m$, that is,\n","\n","$$\n","\\langle s\\rangle = \\mu \\langle v\\rangle = \\mu,\\quad {\\rm var}(s) = \\mu^2{\\rm var}(v) = \\mu^2/m.\n","$$\n","\n","The objective of the speckle filter is to estimate the true signal strength $\\mu$ in a local window which is moved across the whole image. This is a kind of __convolution__, which we will meet later in the form of convolutional neural networks.\n","\n","We'll call the estimated, or filtered, value of the signal strength within the window $\\hat \\mu$ and call the mean observed signal in the window $\\hat s$. \n","\n","Assume there is a linear relation between the contribution to $\\hat \\mu$ from the mean $\\hat s$ and the signal $s$ at the central pixel location,\n","$$\n","\\hat \\mu = a\\hat s + b s.\n","$$\n","We would prefer that in very uniform regions of the SAR scene, $\\hat \\mu \\approx \\hat s$ whereas in rapidly varying regions it is best simply to say $\\hat \\mu \\approx s$.\n","\n","The Lee filter is obtained by minimizing the mean square error \n","\n","$$\n","E = \\langle(\\hat \\mu - \\mu)^2\\rangle = \\langle(a\\hat s + b s -\\mu)^2\\rangle\n","$$\n","\n","with repect to $a$ and $b$. This is done by solving the equations\n","\n","$$\n","{\\partial E\\over \\partial a} = 0, \\quad {\\partial E\\over \\partial b} = 0.\n","$$\n","\n","The result is, see [Canty (2019)](https://www.amazon.com/Analysis-Classification-Change-Detection-Sensing/dp/1138613223/ref=dp_ob_title_bk) p. 193 or [this ESA document](https://earth.esa.int/documents/653194/656796/Speckle_Filtering.pdf):\n","$$\n","\\hat \\mu = \\hat s + {\\hat{\\rm var}(\\mu)\\over \\hat{\\rm var}(s)}(s-\\hat s)\n","$$\n","where $\\hat{\\rm var}(s)$ is the local variance of the observed signal in the window and $\\hat{\\rm var}(\\mu)$ is the local variance of the underlying signal strength. The latter is given by\n","$$\n","\\hat{\\rm var}(\\mu) = {\\hat{\\rm var}(s) -\\hat s^2/m \\over 1+1/m}.\n","$$\n","The $7\\times 7$ window in which mean $\\hat s$ and variance ${\\rm var}(s)$ are evaluated is directionally sensitive and edge-preserving. The masks define the most homogeneous part of the window.\n","\n",""]},{"cell_type":"markdown","metadata":{"id":"cClM-tAWSnuS"},"source":["## Preliminaries"]},{"cell_type":"code","metadata":{"id":"ZdyBUBN4F1IA"},"source":["#@title\n","import ee\n"," \n","# Trigger the authentication flow.\n","ee.Authenticate()\n"," \n","# Initialize the library.\n","ee.Initialize()"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"l1iCQUIDH6wR"},"source":["from google.colab import drive, files\n","drive.mount('/content/drive')"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"X7Z-Enl6Id4n"},"source":["!ls -l /content/drive/MyDrive"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"NqcghIcEIWnZ"},"source":["#@title\n","# Import the Folium library.\n","import folium\n","\n","# Define a method for displaying Earth Engine image tiles to folium map.\n","def add_ee_layer(self, ee_image_object, vis_params, name):\n"," map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)\n"," folium.raster_layers.TileLayer(\n"," tiles = map_id_dict['tile_fetcher'].url_format,\n"," attr = 'Map Data © <a href=\"https://earthengine.google.com/\">Google Earth Engine</a>',\n"," name = name,\n"," overlay = True,\n"," control = True\n"," ).add_to(self)\n","\n","# Add EE drawing method to folium.\n","folium.Map.add_ee_layer = add_ee_layer"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"CTjfQvH1BCgE","cellView":"form"},"source":["#@title Refined Lee filter code\n","# Refined Lee Speckle Filter for S1 images only\n","# Created on 03.03.2020\n","# Transcribed from Guido Lemoine's 2017 JavaScript implementation on the GEE\n","def rl(img):\n","# img must be in natural units, and single band\n","# Set up 3x3 kernels \n"," weights3 = ee.List.repeat(ee.List.repeat(1,3),3)\n"," kernel3 = ee.Kernel.fixed(3,3, weights3, 1, 1, False)\n","\n"," mean3 = img.reduceNeighborhood(ee.Reducer.mean(), kernel3)\n"," variance3 = img.reduceNeighborhood(ee.Reducer.variance(), kernel3)\n","\n","# Use a sample of the 3x3 windows inside a 7x7 window to determine gradients and directions\n"," sample_weights = ee.List([[0,0,0,0,0,0,0], [0,1,0,1,0,1,0],[0,0,0,0,0,0,0], [0,1,0,1,0,1,0], [0,0,0,0,0,0,0], [0,1,0,1,0,1,0],[0,0,0,0,0,0,0]])\n","\n"," sample_kernel = ee.Kernel.fixed(7,7, sample_weights, 3,3, False)\n","\n","# Calculate mean and variance for the sampled windows and store as 9 bands\n"," sample_mean = mean3.neighborhoodToBands(sample_kernel) \n"," sample_var = variance3.neighborhoodToBands(sample_kernel)\n","\n","# Determine the 4 gradients for the sampled windows\n"," gradients = sample_mean.select(1).subtract(sample_mean.select(7)).abs()\n"," gradients = gradients.addBands(sample_mean.select(6).subtract(sample_mean.select(2)).abs())\n"," gradients = gradients.addBands(sample_mean.select(3).subtract(sample_mean.select(5)).abs())\n"," gradients = gradients.addBands(sample_mean.select(0).subtract(sample_mean.select(8)).abs())\n","\n","# And find the maximum gradient amongst gradient bands\n"," max_gradient = gradients.reduce(ee.Reducer.max())\n","\n","# Create a mask for band pixels that are the maximum gradient\n"," gradmask = gradients.eq(max_gradient)\n","\n","# duplicate gradmask bands: each gradient represents 2 directions\n"," gradmask = gradmask.addBands(gradmask)\n","\n","# Determine the 8 directions\n"," directions = sample_mean.select(1).subtract(sample_mean.select(4)).gt(sample_mean.select(4).subtract(sample_mean.select(7))).multiply(1)\n"," directions = directions.addBands(sample_mean.select(6).subtract(sample_mean.select(4)).gt(sample_mean.select(4).subtract(sample_mean.select(2))).multiply(2))\n"," directions = directions.addBands(sample_mean.select(3).subtract(sample_mean.select(4)).gt(sample_mean.select(4).subtract(sample_mean.select(5))).multiply(3))\n"," directions = directions.addBands(sample_mean.select(0).subtract(sample_mean.select(4)).gt(sample_mean.select(4).subtract(sample_mean.select(8))).multiply(4))\n","# The next 4 are the not() of the previous 4\n"," directions = directions.addBands(directions.select(0).Not().multiply(5))\n"," directions = directions.addBands(directions.select(1).Not().multiply(6))\n"," directions = directions.addBands(directions.select(2).Not().multiply(7))\n"," directions = directions.addBands(directions.select(3).Not().multiply(8))\n","\n","# Mask all values that are not 1-8\n"," directions = directions.updateMask(gradmask)\n","\n","# \"collapse\" the stack into a singe band image (due to masking, each pixel has just one value (1-8) in it's directional band, and is otherwise masked)\n"," directions = directions.reduce(ee.Reducer.sum())\n","\n","# var pal = ['ffffff','ff0000','ffff00', '00ff00', '00ffff', '0000ff', 'ff00ff', '000000'];\n","# Map.addLayer(directions.reduce(ee.Reducer.sum()), {min:1, max:8, palette: pal}, 'Directions', false);\n","\n"," sample_stats = sample_var.divide(sample_mean.multiply(sample_mean))\n","\n","# Calculate localNoiseVariance\n"," sigmaV = sample_stats.toArray().arraySort().arraySlice(0,0,5).arrayReduce(ee.Reducer.mean(), [0])\n","\n","# Set up the 7*7 kernels for directional statistics\n"," rect_weights = ee.List.repeat(ee.List.repeat(0,7),3).cat(ee.List.repeat(ee.List.repeat(1,7),4))\n","\n"," diag_weights = ee.List([[1,0,0,0,0,0,0], [1,1,0,0,0,0,0], [1,1,1,0,0,0,0], \n"," [1,1,1,1,0,0,0], [1,1,1,1,1,0,0], [1,1,1,1,1,1,0], [1,1,1,1,1,1,1]]);\n","\n"," rect_kernel = ee.Kernel.fixed(7,7, rect_weights, 3, 3, False)\n"," diag_kernel = ee.Kernel.fixed(7,7, diag_weights, 3, 3, False)\n","\n","# Create stacks for mean and variance using the original kernels. Mask with relevant direction.\n"," dir_mean = img.reduceNeighborhood(ee.Reducer.mean(), rect_kernel).updateMask(directions.eq(1));\n"," dir_var = img.reduceNeighborhood(ee.Reducer.variance(), rect_kernel).updateMask(directions.eq(1))\n","\n"," dir_mean = dir_mean.addBands(img.reduceNeighborhood(ee.Reducer.mean(), diag_kernel).updateMask(directions.eq(2)))\n"," dir_var = dir_var.addBands(img.reduceNeighborhood(ee.Reducer.variance(), diag_kernel).updateMask(directions.eq(2)))\n","\n","# and add the bands for rotated kernels\n"," for i in range(1,4):\n"," dir_mean = dir_mean.addBands(img.reduceNeighborhood(ee.Reducer.mean(), rect_kernel.rotate(i)).updateMask(directions.eq(2*i+1)))\n"," dir_var = dir_var.addBands(img.reduceNeighborhood(ee.Reducer.variance(), rect_kernel.rotate(i)).updateMask(directions.eq(2*i+1)))\n"," dir_mean = dir_mean.addBands(img.reduceNeighborhood(ee.Reducer.mean(), diag_kernel.rotate(i)).updateMask(directions.eq(2*i+2)))\n"," dir_var = dir_var.addBands(img.reduceNeighborhood(ee.Reducer.variance(), diag_kernel.rotate(i)).updateMask(directions.eq(2*i+2)))\n","\n","# \"collapse\" the stack into a single band image (due to masking, each pixel has just one value in it's directional band, and is otherwise masked)\n"," dir_mean = dir_mean.reduce(ee.Reducer.sum())\n"," dir_var = dir_var.reduce(ee.Reducer.sum())\n","\n","# And finally generate the filtered value\n"," varX = dir_var.subtract(dir_mean.multiply(dir_mean).multiply(sigmaV)).divide(sigmaV.add(1.0))\n"," b = varX.divide(dir_var)\n"," result = dir_mean.add(b.multiply(img.subtract(dir_mean))) \n"," return result.arrayFlatten([['sum']])\n","\n","def refinedLee(img):\n"," return ee.Image.cat(rl(img.select(0)),rl(img.select(1)))\n"," \n","if __name__ == '__main__':\n"," pass"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"x3aYgWOGq1pc"},"source":["## SAR time series classification\n","\n","\n"]},{"cell_type":"markdown","metadata":{"id":"0BnDLShSTr7c"},"source":["First of all we grab a time series for an area of interest (aoi) in Saskatchewan, Canada over the 2017 growing season (March to October), see [geojson.io](http://geojson.io/#map=2/20.0/0.0):"]},{"cell_type":"code","metadata":{"id":"rOpuOdJpTpmM"},"source":["geoJSON = {\n"," \"type\": \"FeatureCollection\",\n"," \"features\": [\n"," {\n"," \"type\": \"Feature\",\n"," \"properties\": {},\n"," \"geometry\": {\n"," \"type\": \"Polygon\",\n"," \"coordinates\": [[[-105.10328600000001, 50.24327999999998], \n"," [-104.66649400000001, 50.24327999999998], \n"," [-104.66649400000001, 50.46604134048255], \n"," [-105.10328600000001, 50.46604134048255], \n"," [-105.10328600000001, 50.24327999999998]]]}}]}\n"," \n","coords = geoJSON['features'][0]['geometry']['coordinates']\n","aoi = ee.Geometry.Polygon(coords)\n","\n","collection = ee.ImageCollection('COPERNICUS/S1_GRD_FLOAT') \\\n"," .filterBounds(aoi) \\\n"," .filterDate(ee.Date('2017-03-01'), ee.Date('2017-11-01')) \\\n"," .filter(ee.Filter.eq('transmitterReceiverPolarisation', ['VV','VH'])) \\\n"," .filter(ee.Filter.eq('resolution_meters', 10)) \\\n"," .filter(ee.Filter.eq('instrumentMode', 'IW')) \\\n"," .filter(ee.Filter.eq('relativeOrbitNumber_start',5)) \\\n"," .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING')) \n","collection = collection.sort('system:time_start')\n","\n","collection.size().getInfo()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"WUYkfTQnU4Vl"},"source":["Next, we transform the collection into a multiband image, with and without pre-processing with the refined Lee filter:"]},{"cell_type":"code","metadata":{"id":"JGUMenwTUBOl"},"source":["def get_vvvh(image): \n"," ''' get 'VV' and 'VH' bands from sentinel-1 imageCollection '''\n"," return image.select('VV','VH')\n"," \n","timeseries_rl = collection \\\n"," .map(get_vvvh) \\\n"," .map(refinedLee) \\\n"," .toBands() \\\n"," .clip(aoi) \n","timeseries = collection \\\n"," .map(get_vvvh) \\\n"," .toBands() \\\n"," .clip(aoi) \n","\n","timeseries.bandNames().length().getInfo()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"QglqC-bSVHmM"},"source":["The class lables are conveniently obtained from the GEE archive of the [Canadian AAFC Annual Crop Inventory](https://developers.google.com/earth-engine/datasets/catalog/AAFC_ACI) for the year 2017, and we append them as an additional band (band 39): "]},{"cell_type":"code","metadata":{"id":"IJGTfCPmVBFt"},"source":["crop2017 = ee.ImageCollection('AAFC/ACI') \\\n"," .filter(ee.Filter.date('2017-01-01', '2017-12-01')) \\\n"," .first() \\\n"," .clip(aoi) \\\n"," .float()\n","\n","labeled_timeseries = ee.Image.cat(timeseries, crop2017).float()\n","labeled_timeseries_rl = ee.Image.cat(timeseries_rl, crop2017).float()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"vC6qHGfF3i-D"},"source":["Here is a projection of the class labels and an rgb composite of three VV bands spaced over the time series using Folium. The fact that the colors are intense is encouraging, as it indicates the time series may be discriminating the crop types fairly well:"]},{"cell_type":"code","metadata":{"id":"aP-7xgsrVkeS"},"source":["location = aoi.centroid().coordinates().getInfo()[::-1]\n","\n","timeseries_db = timeseries.log10().multiply(10)\n","timeseries_rl_db = timeseries_rl.log10().multiply(10)\n","\n","rgb = ee.Image.rgb(timeseries_db.select(10),timeseries_db.select(20),timeseries_db.select(30))\n","rgb_rl = ee.Image.rgb(timeseries_rl_db.select(10),timeseries_rl_db.select(20),timeseries_rl_db.select(30))\n","\n","m = folium.Map(location=location, zoom_start=11, height=800, width=1000)\n","\n","m.add_ee_layer(rgb, {'min': [-20], 'max': [0]}, 'rgb')\n","m.add_ee_layer(rgb_rl, {'min': [-20], 'max': [0]}, 'rgb_rl')\n","m.add_ee_layer(labeled_timeseries.select(38), {'min': 0, 'max': 200}, 'crops')\n","\n","m.add_child(folium.LayerControl())"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"arBSr1hZJ93d"},"source":["We can't run TensorFlow on GEE, so to classify the time series with a neural network classifier we will export the filtered and unfiltered labeled time series to GoogleDrive:"]},{"cell_type":"code","metadata":{"id":"bBcrKJbVDqqA"},"source":["drexport = ee.batch.Export.image.toDrive(labeled_timeseries_rl,\n"," description='driveExportTask', \n"," folder = 'gee',\n"," fileNamePrefix='labeled_timeseries_rl',scale=30,maxPixels=1e10)\n","drexport.start()"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"sz_X0y3Ec9Wr"},"source":["drexport1 = ee.batch.Export.image.toDrive(labeled_timeseries,\n"," description='driveExportTask', \n"," folder = 'gee',\n"," fileNamePrefix='labeled_timeseries',scale=30,maxPixels=1e10)\n","drexport1.start()"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"roqBl0gGIori"},"source":["!ls -l /content/drive/MyDrive/gee"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"yOgWc7GtPVG3"},"source":["And then import them into two numpy arrays using _gdal_:"]},{"cell_type":"code","metadata":{"id":"60Ds8T24OTlD"},"source":["from osgeo import gdal\n","from osgeo.gdalconst import GA_ReadOnly, GDT_Byte, GDT_Float32\n","import numpy as np\n","\n","gdal.AllRegister() \n","\n","inDataset = gdal.Open('/content/drive/MyDrive/gee/labeled_timeseries_rl.tif',GA_ReadOnly)\n","cols = inDataset.RasterXSize\n","rows = inDataset.RasterYSize \n","bands = inDataset.RasterCount\n","labeled_timeseries_rl = np.zeros((rows*cols,bands)) \n","for b in range(bands):\n"," band = inDataset.GetRasterBand(b+1)\n"," labeled_timeseries_rl[:,b]=band.ReadAsArray(0,0,cols,rows).ravel() \n","labeled_timeseries_rl = np.nan_to_num(labeled_timeseries_rl) \n","\n","inDataset = gdal.Open('/content/drive/MyDrive/gee/labeled_timeseries.tif',GA_ReadOnly)\n","labeled_timeseries = np.zeros((rows*cols,bands)) \n","for b in range(bands):\n"," band = inDataset.GetRasterBand(b+1)\n"," labeled_timeseries[:,b]=band.ReadAsArray(0,0,cols,rows).ravel() \n","labeled_timeseries = np.nan_to_num(labeled_timeseries) \n","\n","# for later\n","driver = inDataset.GetDriver() \n","inDataset = None\n","m = labeled_timeseries.shape[0]\n","labeled_timeseries.shape"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"0aVeZhhoKsPd"},"source":["The AAFC/ACI thematic maps have 71 different classes. This code generates a dictionary of class names:"]},{"cell_type":"code","metadata":{"id":"FWs8svgLNMec"},"source":["classdict = {'0':'Nc'}\n","filepath = '/content/drive/MyDrive/scripts/AAFC.txt'\n","with open(filepath) as fp:\n"," line = fp.readline()\n"," key = line[:3].replace('\\t','')\n"," value = line[10:].replace('\\t',' ').replace('\\n','')\n"," classdict.update({key:value})\n"," while line:\n"," line = fp.readline()\n"," key = line[:3].replace('\\t','')\n"," value = line[10:].replace('\\t','').replace('\\n','')\n"," classdict.update({key:value})\n","del classdict['']\n","\n","len(classdict)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"aUFxiMI-Pm3l"},"source":["Now we can see which class labels pertain to our region of interest:"]},{"cell_type":"code","metadata":{"id":"YUYiXt56Nl0B"},"source":["classnums = np.unique(labeled_timeseries[:,-1])\n","print(classnums)\n","classnames = str([classdict[str(int(cn))] for cn in classnums])\n","classnames"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"DQfZVMRUPvim"},"source":["In order to train the neural network we have to renumber the labels consecutively from 0:"]},{"cell_type":"code","metadata":{"id":"x4BtyCekP1nR"},"source":["i=0\n","labels = labeled_timeseries[:,-1]\n","for c in classnums:\n"," labels = np.where(labels==c,i,labels)\n"," i += 1 \n","labels = np.array(labels,dtype=np.uint8) \n","n_classes = len(np.unique(labels))\n","print(np.unique(labels))"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"BEhtLLnnUM9-"},"source":["Here we write the class labels image to the Drive"]},{"cell_type":"code","metadata":{"id":"ILno5IcRQS1O"},"source":["outDataset = driver.Create('/content/drive/MyDrive/gee/labels.tif',cols,rows,1,GDT_Byte)\n","outBand = outDataset.GetRasterBand(1)\n","outBand.WriteArray(np.reshape(labels,(rows,cols)))\n","outBand.FlushCache()\n","outDataset = None"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"WkvhJvTX2BcW"},"source":[" and display it using the script dispms.py (get help with the -h option):"]},{"cell_type":"code","metadata":{"id":"rrRXnTx72Ktd"},"source":["%run /content/drive/MyDrive/scripts/dispms -f /content/drive/MyDrive/gee/labels.tif -c "],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"I3fIZmh1VBRm"},"source":["Now we simulate ground truth by taking a (generous) random subset of 50k training pixels, one from the unfiltered and one from the filtered time series:"]},{"cell_type":"code","metadata":{"id":"X7lS_WSyVM0x"},"source":["# random subset for training\n","np.random.seed(321)\n","n = 50000\n","idx = np.random.permutation(m)[0:n]\n","# training vectors (x100)\n","Xs = labeled_timeseries[idx,:-1]*100\n","Xs_rl = labeled_timeseries_rl[idx,:-1]*100\n","\n","# one hot encoded class labels\n","Ls = np.array(labels[idx],dtype=np.int)\n","ls = tf.one_hot(Ls,n_classes)\n","print(ls[0:5,:]) "],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"keOciPVBVaVT"},"source":["Now we instantiate the neural network class Ffn() with two hidden layers, each with 40 neurons:"]},{"cell_type":"code","metadata":{"id":"7BrBUw1LViB7"},"source":["classifier = Ffn(Ls=[40,40], K=n_classes, learning_rate=0.002)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"FMCET02E7QIH"},"source":["And train it on the time series:"]},{"cell_type":"code","metadata":{"id":"dMJoqbc37OLH"},"source":["%%time\n","classifier.train(Xs, ls, epochs=200)"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"rqXXdGG4Yd54"},"source":["classifier.history()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"QUBy9VCChFsz"},"source":["So the accuracy is around 82%. Now let's try with the filtered image series, without re-setting the classifier (a foretaste of _transfer learning_, the subject of the next part of the course):"]},{"cell_type":"code","metadata":{"id":"BAYIeJJdhXJQ"},"source":["%%time\n","classifier.train(Xs_rl, ls, epochs=100)"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"2TwNiGnaiX4Z"},"source":["classifier.history()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"tvrtHbioaIZw"},"source":["Test the classifier with all of the training data (we should really have reserved a test set for this):"]},{"cell_type":"code","metadata":{"id":"e71GQ4BxaNAM"},"source":["classifier.test(Xs_rl,ls)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"oERdhzSuYy7e"},"source":["This is a considerable improvement, so we'll use this result to classify (predict) the entire image:"]},{"cell_type":"code","metadata":{"id":"BRs6Axe0ZFPu"},"source":["cls, probs = classifier.classify(labeled_timeseries_rl[:,0:-1]*100) \n","# for later display:\n","cls[0]=1\n","cls[1]=n_classes-1\n","\n","probs.shape"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"Ve_oVjSCZr9D"},"source":["Now write the thematic map and the class probabilities images to Drive:"]},{"cell_type":"code","metadata":{"id":"FiqXj1mbZx0g"},"source":["# write the class image to disk\n","outDataset = driver.Create('/content/drive/MyDrive/gee/timeseries_class.tif',cols,rows,1,GDT_Byte)\n","outBand = outDataset.GetRasterBand(1)\n","outBand.WriteArray(np.reshape(cls,(rows,cols)))\n","outBand.FlushCache()\n","outDataset = None\n","# write the class probabilities to disk\n","bands = probs.shape[1]\n","probs = np.byte(probs*255)\n","outDataset = driver.Create('/content/drive/MyDrive/gee/timeseries_probs.tif',cols,rows,bands,GDT_Float32)\n","for b in range(bands):\n"," outBand = outDataset.GetRasterBand(b+1)\n"," outBand.WriteArray(np.reshape(probs[:,b],(rows,cols)))\n"," outBand.FlushCache()\n","outDataset = None"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"XApieGNoaXvI"},"source":["Compare the classified image with the AAFC/ACI thematic map:"]},{"cell_type":"code","metadata":{"id":"EWL8GMq5ac5D"},"source":["%run /content/drive/MyDrive/scripts/dispms -f /content/drive/MyDrive/gee/timeseries_class.tif -c -F /content/drive/MyDrive/gee/labels.tif -C "],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"0Sf6aFPea7n7"},"source":["Looks good, but we can \"pretty it up\" some more with "]},{"cell_type":"markdown","metadata":{"id":"6Kxbxo5zsMCC"},"source":["## Probabilistic Label Relaxation:\n","\n","With a heuristic ansatz [Richards (2012)](https://www.amazon.de/Remote-Sensing-Digital-Image-Analysis-ebook/dp/B00A9YH460/ref=sr_1_3?__mk_de_DE=%C3%85M%C3%85%C5%BD%C3%95%C3%91&dchild=1&keywords=richards+remote+sensing&qid=1611397419&sr=8-3) the output probability $m_k$ for an input pixel vector $g$ can be corrected to\n","$$\n","m_k \\to {m_k Q_k\\over \\sum_j m_jQ_j},\n","$$\n","where $Q$ is a neighbourhood function which can be determined by the contextual information in the classified image, see [Canty (2019)]((https://www.amazon.com/Analysis-Classification-Change-Detection-Sensing/dp/1138613223/ref=dp_ob_title_bk), p. 294."]},{"cell_type":"code","metadata":{"id":"g7rm5yhpojDs"},"source":["%run /content/drive/MyDrive/scripts/dispms -f /content/drive/MyDrive/gee/timeseries_probs.tif -p [12,14,16] -e 3"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"_KvSUH56bSGQ"},"source":["%run /content/drive/MyDrive/scripts/plr /content/drive/MyDrive/gee/timeseries_probs.tif"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"0fAck2vwcyWL"},"source":["%run /content/drive/MyDrive/scripts/dispms -f /content/drive/MyDrive/gee/timeseries_probs_plr.tif -c -F /content/drive/MyDrive/gee/timeseries_class.tif -C"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"XpgTT-9yddHw"},"source":["### Project\n","Try to apply a similar procedure to a Sentinel-2 time series."]},{"cell_type":"markdown","metadata":{"id":"GUR6IUKR-lv4"},"source":["### Outlook\n","In the next part of the course we'll leave the subject of polarimetric SAR time series and have a look at \"deep learning\" classification of Sentinel-2 imagery with convolutional neural networks."]}]} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Looking forward to the next one