<jupyter_start><jupyter_text>Google Earth Engine, Machine Learning  and Multitemporal Sentinel-1 Processing  ZFL Short Course, April 2020 Part 1. The GEE "Earth Engine is optimized for data processing and transformation, by co-locating the data with the CPUs.  It's good at turning lots of big inputs into a much smaller set of results; one canonical example is generating temporal composites of band indices, where you might be processing over an entire year's worth of data to produce a single best-pixel image of some transformation, and then maybe producing a land-cover classification from that and computing areal statistics per land cover class.The whole point of Earth Engine is to keep users from having to manage downloading a bunch of data and pending time dealing with storage, file formats, network bandwidth, CPUs, package installation, parallelism and batch queues as much as possible.  So it is not optimized for downloading lots of data, and in fact, you're likely to find that getting data out is likely to be the biggest bottleneck.  It's not that you can't do it, it's just that that's not what it's best at."-<NAME>  Four ways to access and program the GEE 1. Work entirely in the code editor with the JavaScript APIhttps://code.earthengine.google.com/   2. Work with the Python API on the Web (Colab)https://colab.research.google.com/notebooks/welcome.ipynb 3. Work with the Python API locally (conda)https://developers.google.com/earth-engine/python_install-conda.html 4. Use my Docker container for this course (Jupyter Notebooks)https://github.com/mortcanty/EESARDocker  Why not use Colab? Answer: Can't import my auxil package (without putting it in PyPI).  Can't run external Python scripts Why not use conda? You can if you like. Good luck! Why not use Jupyterlab? Answer: Jupyter widgets, ipyleaflet  The Python API on a Jupyter Notebook Served from a Docker Container docker run -d -p 8888:8888 --name=eesar mort/eesardocker Point your browser to localhost:8888 Open a terminal window and authenticate  The command line interface<jupyter_code>!earthengine<jupyter_output><empty_output><jupyter_text>Enable inline graphics<jupyter_code>%matplotlib inline<jupyter_output><empty_output><jupyter_text>Accessing the data catalogue<jupyter_code>import ee, time

ee.Initialize()

try:
    point = ee.Geometry.Point([8.5,50.0]) 
    start = ee.Date('2016-05-01')
    finish = ee.Date('2016-08-01')    
    collection = ee.ImageCollection('COPERNICUS/S1_GRD') \
                .filterBounds(point) \
                .filterDate(start, finish) \
                .filter(ee.Filter.eq('transmitterReceiverPolarisation', ['VV','VH'])) \
                .filter(ee.Filter.eq('resolution_meters', 10)) \
                .filter(ee.Filter.eq('instrumentMode', 'IW')) \
                .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
    count = collection.size().getInfo()
    if count==0:
        raise ValueError('No images found')   
    image = collection.first()
    timestamp = ee.Date(image.get('system:time_start')).getInfo()
    timestamp = time.gmtime(int(timestamp['value'])/1000)
    timestamp = time.strftime('%c', timestamp) 
    systemid = image.get('system:id').getInfo()
    projection = image.select(0).projection().getInfo()['crs']    
    print('system id: %s'%systemid)
    print('acquisition time: %s'%timestamp)
    print('projection: %s'%projection)
except Exception as e:
    print('An error occurred: %s'%e)<jupyter_output>system id: COPERNICUS/S1_GRD/S1A_IW_GRDH_1SDV_20160504T171539_20160504T171604_011112_010BED_80AF
acquisition time: Wed May  4 17:15:39 2016
projection: EPSG:32632<jupyter_text>Exporting data<jupyter_code>maxLat = 50.06526459341422
minLon = 8.477334975832491
minLat = 50.01013697421883
maxLon = 8.633890151613743
rect = ee.Geometry.Rectangle(minLon,minLat,maxLon,maxLat)
exportname = 'users/mortcanty/'+systemid.replace('/','-')
assexport = ee.batch.Export.image.toAsset(image.clip(rect).select('B2','B3','B4','B8'),
                                          description='assetExportTask', 
                                          assetId=exportname,scale=10,maxPixels=1e9)
assexportid = str(assexport.id)
print('****Exporting to Assets, task id: %s '%assexportid)
assexport.start()<jupyter_output><empty_output><jupyter_text>Programming the GEE  Mapping<jupyter_code>alist = ee.List([1,2,3,4])

def squareit(x):
    x = ee.Number(x)
    return x.multiply(x)

squaredlist = alist.map(squareit)
print(squaredlist.getInfo())<jupyter_output><empty_output><jupyter_text>Iteration<jupyter_code># Pure Python
sum = 0
for i in range(10):
    sum += i
print(sum)
# numpy
import numpy as np

print(np.sum(range(10)))
# ee Python API
def sum(current,prev):
    prev = ee.Number(prev)
    return prev.add(current)

seq = ee.List.sequence(0,9)
first = ee.Number(0)
result = seq.iterate(sum,first)
print(result.getInfo())<jupyter_output>45.0<jupyter_text>Reducers<jupyter_code>seq.reduce(ee.Reducer.sum()).getInfo()
# Image.reduceRegion example
#
# Computes a simple reduction over a region of an image.  A reduction
# is any process that takes an arbitrary number of inputs (such as
# all the pixels of an image in a given region) and computes one or
# more fixed outputs.  The result is a dictionary that contains the
# computed values, which in this example is the maximum pixel value
# in the region.
#
# This example shows how to print the resulting dictionary to the
# console, which is useful when developing and debugging your
# scripts, but in a larger workflow you might instead use the
# Dicitionary.get() function to extract the values print np.max(band1)you need from the
# dictionary for use as inputs to other functions.

# The input image to reduce, in this case an SRTM elevation map.
image = ee.Image('srtm90_v4')

# The region to reduce within.
poly = ee.Geometry.Rectangle(-109.05, 41, -102.05, 37)

# Reduce the image within the given region, using a reducer that
# computes the max pixel value.  We also specify the spatial
# resolution at which to perform the computation, in this case 200
# meters.
max = image.reduceRegion(ee.Reducer.max(), poly, 200)

# Print the result to the console.
print(max.getInfo())<jupyter_output>{'elevation': 4371}<jupyter_text>Ipyleaflet In Python, we don't have the benefit of the GEE Payground (code editor) so we need an ersatz<jupyter_code>from ipyleaflet import Map, GeoJSON, TileLayer
import json
import os
import requests
import ee

ee.Initialize()

def GetTileLayerUrl(ee_image_object):
    map_id = ee.Image(ee_image_object).getMapId()
    tile_url_template = "https://earthengine.googleapis.com/map/{mapid}/{{z}}/{{x}}/{{y}}?token={token}"
    return tile_url_template.format(**map_id)

if not os.path.exists('europe_110.geo.json'):
  url = 'https://github.com/jupyter-widgets/ipyleaflet/raw/master/examples/europe_110.geo.json'
  r = requests.get(url)
  with open('europe_110.geo.json', 'w') as f:
    f.write(r.content.decode("utf-8"))

m = Map(center=(50.6252978589571, 0.34580993652344), zoom=3)
m
with open('europe_110.geo.json', 'r') as f:
  data = json.load(f)
geo_json = GeoJSON(data=data, style = {'color': 'green', 'opacity':1, 'weight':1.9, 'dashArray':'9', 'fillOpacity':0.1})
m.add_layer(geo_json)
rect = ee.Geometry.Rectangle( [5.9985351562500009, 50.938486572440667, 6.0946655273437509, 50.973953836311068])
avimg = ee.ImageCollection('COPERNICUS/S1_GRD') \
                        .filterBounds(rect) \
                        .filterDate(ee.Date('2018-04-01'), ee.Date('2018-09-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('orbitProperties_pass', 'ASCENDING')) \
                        .sort('system:time_start') \
                        .mean() \
                        .select(0) \
                        .visualize(min=-15,max=1)
m.add_layer(TileLayer(url=GetTileLayerUrl(avimg)))<jupyter_output><empty_output><jupyter_text>Part 2. Polarimetric SAR   Vector and matrix representationsA fully polarimetric SAR measures a$2\times 2$ _scattering matrix_ $S$  at each resolution cell on the ground.The scattering matrix relates the incident and the backscatteredelectric fields $E^i$ and $E^b$ according to$$\pmatrix{E_h^b \cr E_v^b}=\pmatrix{S_{hh} & S_{hv}\cr S_{vh} & S_{vv}}\pmatrix{E_h^i \cr E_v^i}.$$Here $E_h^{i(b)}$ and $E_v^{i(b)}$ denote the horizontal and vertical components of the incident (backscattered)oscillating electric fields directly at the target. These can be deduced from the transmitted and receivedradar signals via the so-called _far field_ approximations.If both horizontally and vertically polarized radar pulses areemitted and discriminated then they determine, from the above Equation, the four complex scattering matrix elements.The per-pixel polarimetric information in the scattering matrix $S$, under the assumptionof reciprocity ($S_{hv} = S_{vh}$), can then be expressed as a three-component complex vector$$s = \pmatrix{S_{hh}\cr \sqrt{2}S_{hv}\cr S_{vv}},$$where the $\sqrt{2}$ ensures that the total intensity (received signal power) is consistent. The total intensity is referred to as the _span_ and is the complex inner product of the vector $s$,$${\rm span} = s^\top s = |S_{hh}|^2 + 2|S_{hv}|^2 + |S_{vv}|^2.$$This is a real number and the corresponding gray-scale image is called the _span image_. The observation vector $s$ can be shown to be a realization of a complex multivariate normal random variable. An equivalent and often preferred representation is in terms of the *coherency vector* or *Pauli decomposition*$$k = {1\over\sqrt{2}}\pmatrix{S_{hh} + S_{vv}\cr S_{hh} - S_{vv} \cr 2S_{hv}}.$$The polarimetric signal is  can also be represented by taking the complex outer product of $s$ with itself:$$C = s s^\top = \pmatrix{ |S_{hh}|^2 & \sqrt{2}S_{hh}S_{hv}^* & S_{hh}S_{vv}^* \cr                                     \sqrt{2}S_{hv}S_{hh}^* & 2|S_{hv}|^2 & \sqrt{2}S_{hv}S_{vv}^* \cr                                     S_{vv}S_{hh}^* & \sqrt{2}S_{vv}S_{hv}^* & |S_{vv}|^2 }.$$The diagonal elements of $C$ are real numbers, with span $= {\rm tr}(C)$, and the off-diagonalelements are complex. This matrix representation contains all of the information in the polarized signal. Multi-lookingThe matrix $C$ can be averaged over the number of looks (number of adjacent cells used to average out the effect of speckle) to give an estimate of the __covariance matrix__ of each multi-look pixel:$$\bar{C}  ={1\over m}\sum_{\nu=1}^m  s(\nu) s(\nu)^\top = \langle  s s^\top \rangle = \pmatrix{ \langle |S_{hh}|^2\rangle & \langle\sqrt{2}S_{hh}S_{hv}^*\rangle & \langle S_{hh}S_{vv}^*\rangle \cr\langle\sqrt{2} S_{hv}S_{hh}^*\rangle & \langle 2|S_{hv}|^2\rangle & \langle\sqrt{2}S_{hv}S_{vv}^*\rangle \cr\langle S_{vv}S_{hh}^*\rangle & \langle\sqrt{2}S_{vv}S_{hv}^*\rangle & \langle |S_{vv}|^2\rangle },$$where $m$ is the number of looks. This matrix (or alternatively the equivalent __coherency matrix__ $\langle  k k^\top \rangle$) can be generated with the Sentinel-1 Toolbox.  Rewriting the first of the above equalities,$$m\bar{C} = \sum_{\nu=1}^m  s(\nu) s(\nu)^\top =:  x.$$This quantity $x$ is a realization of a __complex random matrix__ which is known to havea complex Wishart distribution with $N\times N$ covariance matrix $\Sigma$ (here $N=3$) and $m$ degrees of freedom:$$p_{W_c}( x) ={|x|^{(m-N)}\exp(-{\rm tr}(\Sigma^{-1} x)) \over\pi^{N(N-1)/2}|\Sigma|^{m}\prod_{i=1}^N\Gamma(m+1-i)},\quad m \ge N,$$provided that the vectors $s(\nu)$ are independent and drawn from a complex multivariate normal distribution. Dual and single polarimetric imageryThe scattering vector given above corresponds to so-called full, or _quad polarimetric_ SAR.Satellite-based SAR sensors often operate in reduced, power-saving polarization modes, emitting only one polarization and receivingtwo (dual polarization) or one (single polarization). The look-averaged covariance matrices are reduced in dimensioncorrespondingly. For instance for dual polarization with horizontal transmission and horizontal and vertical reception,$$\bar{C} = \pmatrix{ \langle |S_{hh}|^2\rangle & \langle S_{hh}S_{hv}^*\rangle \cr\langle S_{hv}S_{hh}^*\rangle & \langle |S_{hv}|^2\rangle },$$and, for single polarization and horizontal transmission/reception, we get simply the intensity image$$\bar{I} = \langle |S_{hh}|^2\rangle \quad {\rm or} \quad \bar{I} = \langle |S_{vv}|^2\rangle.$$A special case, relevant to Sentinel-1 data on the GEE is vertical transmission, vertical and horizontal reception without including the off diagonal complex element:$$\bar{C} = \pmatrix{ \langle |S_{vv}|^2\rangle & 0 \cr 0 & \langle |S_{vh}|^2\rangle }$$referred to as dualpol, diagonal only.  Multi-temporal data: the omnibus testThe following is discussion is based on Conradsen et al (2016).As we have seen, we can represent a pixel vector in an $m$ look-averagedpolSAR image in covariance matrix format by $\bar C$, where$$m\bar C =  x = \sum_{\nu=1}^m  s(\nu) s(\nu)^\top$$is a realization of a random matrix $X$ with a complex Wishart distribution.In order to derive a change detection procedure for two polarimetric SAR images, we formulate a statistical test. For each pixel in $k$ images, define the null (or no-change) simple  hypothesis as$$H_0:\quad \Sigma_1 = \Sigma_2 = \dots = \Sigma_k = \Sigma,$$and the alternative composite hypothesis as$$H_1:\quad \Sigma_i \ne \Sigma_j $$for at least one pair $i,j$Then the __likelihood ratio test__ has the critical region for rejection of the no-change hypothesis$$Q = k^{kNm}{ |x_1|^m |x_2|^m\cdots |x_k|^m \over |x_1 + x_2 + \dots x_k|^{km} } \le t$$where $t$ is some small number and $N$ is the dimension of $x$.Taking logarithms:$$\ln Q = m(Nk\ln k + \sum_{i=1}^k \ln|x_i| - k \ln|x|) \le t'$$where $x = \sum_{i=1}^k x_i$.Finally, we have the following approximation for the statistical distribution of the test statistic $Q$:$${\rm prob}(-2\log Q\le z) \simeq P_{\chi^2;(k-1)f}(z)$$where $f$ is the number of parameters requred to specify $x$: 9 for quadpol, 4 for dualpol and 2 for dualpol diagonal only.In practice we choose a significance level $\alpha$, e.g., $\alpha = 0.01$, and define the p-value as$$p = 1-{\rm prob}(-2\log Q\le z)$$and interpret all pixels with $p<\alpha$ as change.   Sequential omnibus testFurthermore this test can be factored into a sequence of  tests involving hypotheses of the form$\Sigma_1 = \Sigma_2$ against $\Sigma_1 \ne \Sigma_2$,$\Sigma_1 = \Sigma_2 = \Sigma_3$ against $\Sigma_1 = \Sigma_2 \ne \Sigma_3$, and so forth. For example, to test$$H_{0j}: \quad \Sigma_1 = \Sigma_2 = \dots \Sigma_{j-1} = \Sigma_j$$against$$H_{0j}: \quad \Sigma_1 = \Sigma_2 = \dots \Sigma_{j-1} \ne \Sigma_j$$the likelihood ratio test statstic is $R^1_j$ given by$$\ln R^1_j = m\big[N(j\ln j - (j-1)\ln (j-1)) + (j-1)\ln \big|\sum_{i=1}^{j-1} x_i\big| + \ln|x_j| - j\ln \big|\sum_{i=1}^{j} x_i\big|\ \big] \quad j=2\dots k.$$and$${\rm prob}(-2\log R^1_j\le z) \simeq P_{\chi^2;f}(z)$$Now suppose that we conclude $\Sigma_1\ne \Sigma_2$. Then we can continue to look for additional changes by restarting the tests at $j=3$,$$\ln R^2_j = m\big[N(j\ln j - (j-1)\ln (j-1)) + (j-1)\ln \big|\sum_{i=2}^{j-1} x_i\big| + \ln|x_j| - j\ln \big|\sum_{i=2}^{j} x_i\big|\ \big] \quad j=3\dots k.$$For a series of, say, 5 images we therefore have, __for each pixel__, the following tests to consider$$\matrix{R^1_2 & R^1_3 & R^1_4 & R^1_5 \cr      & R^2_3 & R^2_4 & R^2_5 \cr      &       & R^3_4 & R^3_5 \cr      &       &       & R^4_5 }$$<jupyter_code>rect = ee.Geometry.Rectangle( [5.9985351562500009, 50.938486572440667, 6.0946655273437509, 50.973953836311068])
collection = ee.ImageCollection('COPERNICUS/S1_GRD') \
                        .filterBounds(rect) \
                        .filterDate(ee.Date('2018-04-01'), ee.Date('2018-10-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('orbitProperties_pass', 'ASCENDING')) \
                        .sort('system:time_start')                                     
acquisition_times = ee.List(collection.aggregate_array('system:time_start')).getInfo() 
count = len(acquisition_times)
print(count)
from auxil import eeWishart
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()

def clipList(current,prev):
    ''' clip a list of images '''
    imlist = ee.List(ee.Dictionary(prev).get('imlist'))
    rect = ee.Dictionary(prev).get('rect')    
    imlist = imlist.add(ee.Image(current).clip(rect))
    return ee.Dictionary({'imlist':imlist,'rect':rect})

pcollection = collection.map(get_vvvh)

# get the list of images and clip to roi
pList = pcollection.toList(count)   
first = ee.Dictionary({'imlist':ee.List([]),'rect':rect}) 
imList = ee.Dictionary(pList.iterate(clipList,first)).get('imlist')  

# run the algorithm ------------------------------------------   
result = ee.Dictionary(eeWishart.omnibus(imList,0.0001,4.4,True))
# ------------------------------------------------------------      

cmap = ee.Image(result.get('cmap')).byte()   
smap = ee.Image(result.get('smap')).byte()
fmap = ee.Image(result.get('fmap')).byte()  
bmap = ee.Image(result.get('bmap')).byte()
import IPython.display as disp
jet = 'black,blue,cyan,yellow,red'
url = fmap.getThumbURL({'min':0,'max':15,'palette':jet})
disp.Image(url=url,width = 1000)
from auxil.eeSar_seq import run
run()<jupyter_output><empty_output><jupyter_text>Dual Polarimtric Imagery from the Copernicus Open Access Hub<jupyter_code>from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt
import getpass
passwd = getpass.getpass()
api = SentinelAPI('mortcanty',passwd,'https://scihub.copernicus.eu/dhus')
# http://geojson.io/#map=12/50.8866/6.5176
footprint = geojson_to_wkt({
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              6.3775634765625,
              50.81461241998999
            ],
            [
              6.65771484375,
              50.81461241998999
            ],
            [
              6.65771484375,
              50.958426723359935
            ],
            [
              6.3775634765625,
              50.958426723359935
            ],
            [
              6.3775634765625,
              50.81461241998999
            ]
          ]
        ]
      }
    }
  ]
})
products = api.query(footprint,
                     date=('20190501','20190507'),
                     platformname='Sentinel-1',
                     polarisationmode='VV VH',
                     sensoroperationalmode='IW',
                     orbitdirection='Ascending',
                     relativeorbitnumber='15',
                     producttype='SLC')
for key, value in products.items(): 
    print(key,value['title'])
for key, value in products.items(): 
    product_info = api.get_product_odata(key)
    if product_info['Online']:
        print('Product %s is online. Starting download.'%key)
        api.download(key,directory_path='imagery')
    else:
        print('Product %s is not online.'%key)
api.download_all(products,directory_path='/home/mort/python/EESARDocker/src/imagery')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13826512590200332, 'fraction_numerical': 0.03981430085280315, 'mean_word_length': 3.5527222605099933, 'pattern_counts': {'":': 7, '<': 40, '<?xml version=': 0, '>': 39, 'https://': 7, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '292714', 'license': 'MIT', 'n_tokens_mistral': 6982, 'n_tokens_neox': 6317, 'n_words': 1843, 'path': 'src/zfl2020_1.ipynb', 'repo_name': 'jrmanceraf1421/sar', 'size': '19824'}