Calculate a weighted centered image $(y)_{ij}$ from image $(x)_{ij},\ \ i=1\dots N,\ \ j=1\dots n$, with weights $w_j,\ \ j=1\dots n$, where $N$ is the number of bands and $n$ is the number of pixels:
$$ \bar x_i = {1 \over \sum_{j=1}^n w_j} \sum_{j=1}^n w_j x_{ij},\quad i = 1\dots N $$$$ y_{ij} = x_{ij}-\bar x_i,\quad i = 1\dots N,\ j=1\dots n $$Calculate the weighted covariance matrix $(c)_{k\ell},\ \ k,\ell = 1\dots N$, of a weighted centered image $(y)_{ij}$:
$$ c_{k,\ell} = {1\over \sum_{i=1}^n w_i}\sum_{j=1}^n w_j y_{kj}y_{\ell j} = {1\over \sum_{i=1}^n w_i}\sum_{j=1}^n \sqrt{w_j} y_{kj}\sqrt{w_j}y_{\ell j} $$
In [1]:
import ee
ee.Initialize()
In [2]:
def covarw(image, weights, scale=30, maxPixels=1e9):
'''Return the weighted centered image and its weighted covariance matrix'''
bandNames = image.bandNames()
N = bandNames.length()
weightsImage = image.multiply(ee.Image.constant(0)).add(weights)
means = image.addBands(weightsImage) \
.reduceRegion(ee.Reducer.mean().repeat(N).splitWeights(),scale=scale,maxPixels=maxPixels) \
.toArray() \
.project([1])
centered = image.toArray().subtract(means)
B1 = centered.bandNames().get(0)
b1 = weights.bandNames().get(0)
nPixels = ee.Number(centered.reduceRegion(ee.Reducer.count(), scale=scale, maxPixels=maxPixels).get(B1))
sumWeights = ee.Number(weights.reduceRegion(ee.Reducer.sum(), scale=scale, maxPixels=maxPixels).get(b1))
covw = centered.multiply(weights.sqrt()) \
.toArray() \
.reduceRegion(ee.Reducer.centeredCovariance(), scale=scale, maxPixels=maxPixels) \
.get('array')
covw = ee.Array(covw).multiply(nPixels).divide(sumWeights)
return (centered.arrayFlatten([bandNames]), covw)
In [4]:
import numpy as np
minlon = -116.117
minlat = 36.964
maxlon = -115.920
maxlat = 37.109
rect = ee.Geometry.Rectangle(minlon,minlat,maxlon,maxlat)
image = ee.Image('LT5_L1T/LT50400341985097XXX04') \
.select('B1','B2','B3','B4') \
.clip(rect)
npixels = ee.Number(image.select(0).reduceRegion(ee.Reducer.count(), scale=30, maxPixels=1e9).get('B1'))
print 'number of pixels'
print npixels.getInfo()
# equal weights
weights = image.select(0).multiply(0.0).add(1.0)
_, covw = covarw(image,weights)
print 'unweighted covariance matrix, covarw()'
print np.array(ee.Array(covw).getInfo())
# should be same as ordinary covariance
cov = image.toArray().reduceRegion(ee.Reducer.covariance(),scale=30,maxPixels=1e9)
print 'unweighted covariance matrix, ee.Reducer.covariance()'
print np.array(cov.get('array').getInfo())
# different weights (just the pixel values themselves)
weights = image.float().select(0)
_, covw = covarw(image,weights)
print 'weighted covariance matrix, covarw()'
print np.array(ee.Array(covw).getInfo())
In [17]:
gdexport = ee.batch.Export.image.toDrive(image,
description='driveExportTask',
folder = 'EarthEngineImages',
fileNamePrefix='image',scale=30,maxPixels=1e9)
gdexport.start()
gdexport = ee.batch.Export.image.toDrive(weights,
description='driveExportTask',
folder = 'EarthEngineImages',
fileNamePrefix='weights',scale=30,maxPixels=1e9)
gdexport.start()
In [5]:
import gdal,sys
from osgeo.gdalconst import GA_ReadOnly
inDataseti = gdal.Open('image.tif',GA_ReadOnly)
cols = inDataseti.RasterXSize
rows = inDataseti.RasterYSize
bands = inDataseti.RasterCount
tmp = inDataseti.GetRasterBand(1).ReadAsArray(0,0,cols,rows).astype(float).ravel()
idx = np.where(tmp>0)
npixels = np.size(idx)
print 'number of pixels'
print npixels
print 'unweighted covariance matrix'
G = np.zeros((npixels,bands))
for b in range(bands):
band = inDataseti.GetRasterBand(b+1)
tmp = band.ReadAsArray(0,0,cols,rows).astype(float).ravel()
tmp = tmp[idx]
G[:,b] = tmp
C = np.cov(G,rowvar=0)
print C
inDatasetw = gdal.Open('weights.tif',GA_ReadOnly)
cols = inDatasetw.RasterXSize
rows = inDatasetw.RasterYSize
weights = inDatasetw.GetRasterBand(1).ReadAsArray(0,0,cols,rows).astype(float).ravel()
weights = weights[idx]
print 'weighted covariance matrix'
C = np.cov(G,rowvar=0,aweights=weights)
print C
In [ ]: