Skip to content

covariates module

indices

Source code in src\rlcms\covariates.py 
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
class indices():

	def __init__(self):



		# list with functions to call for each index
		self.functionList = {"ND_blue_green" : self.ND_blue_green, \
							 "ND_blue_red" : self.ND_blue_red, \
							 "ND_blue_nir" : self.ND_blue_nir, \
							 "ND_blue_swir1" : self.ND_blue_swir1, \
							 "ND_blue_swir2" : self.ND_blue_swir2, \
							 "ND_green_red" : self.ND_green_red, \
							 "ND_green_nir" : self.ND_green_nir, \
							 "ND_green_swir1" : self.ND_green_swir1, \
							 "ND_green_swir2" : self.ND_green_swir2, \
							 "ND_red_swir1" : self.ND_red_swir1, \
							 "ND_red_swir2" : self.ND_red_swir2, \
							 "ND_nir_red" : self.ND_nir_red, \
							 "ND_nir_swir1" : self.ND_nir_swir1, \
							 "ND_nir_swir2" : self.ND_nir_swir2, \
							 "ND_swir1_swir2" : self.ND_swir1_swir2, \
							 "R_swir1_nir" : self.R_swir1_nir, \
							 "R_red_swir1" : self.R_red_swir1, \
							 "EVI" : self.EVI, \
							 "SAVI" : self.SAVI, \
							 "IBI" : self.IBI}


	def addAllTasselCapIndices(self,img): 
		""" Function to get all tasselCap indices """

		def getTasseledCap(img):
			"""Function to compute the Tasseled Cap transformation and return an image"""

			coefficients = ee.Array([
				[0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863],
				[-0.2848, -0.2435, -0.5436, 0.7243, 0.0840, -0.1800],
				[0.1509, 0.1973, 0.3279, 0.3406, -0.7112, -0.4572],
				[-0.8242, 0.0849, 0.4392, -0.0580, 0.2012, -0.2768],
				[-0.3280, 0.0549, 0.1075, 0.1855, -0.4357, 0.8085],
				[0.1084, -0.9022, 0.4120, 0.0573, -0.0251, 0.0238]
			])

			bands=ee.List(['blue','green','red','nir','swir1','swir2'])

			# Make an Array Image, with a 1-D Array per pixel.
			arrayImage1D = img.select(bands).toArray()

			# Make an Array Image with a 2-D Array per pixel, 6x1.
			arrayImage2D = arrayImage1D.toArray(1)

			componentsImage = ee.Image(coefficients).matrixMultiply(arrayImage2D).arrayProject([0]).arrayFlatten([['brightness', 'greenness', 'wetness', 'fourth', 'fifth', 'sixth']]).float()

			# Get a multi-band image with TC-named bands.
			return img.addBands(componentsImage);	


		def addTCAngles(img):

			""" Function to add Tasseled Cap angles and distances to an image. Assumes image has bands: 'brightness', 'greenness', and 'wetness'."""

			# Select brightness, greenness, and wetness bands	
			brightness = img.select('brightness')
			greenness = img.select('greenness')
			wetness = img.select('wetness')

			# Calculate Tasseled Cap angles and distances
			tcAngleBG = brightness.atan2(greenness).divide(math.pi).rename(['tcAngleBG'])
			tcAngleGW = greenness.atan2(wetness).divide(math.pi).rename(['tcAngleGW'])
			tcAngleBW = brightness.atan2(wetness).divide(math.pi).rename(['tcAngleBW'])
			tcDistBG = brightness.hypot(greenness).rename(['tcDistBG'])
			tcDistGW = greenness.hypot(wetness).rename(['tcDistGW'])
			tcDistBW = brightness.hypot(wetness).rename(['tcDistBW'])
			img = img.addBands(tcAngleBG).addBands(tcAngleGW).addBands(tcAngleBW).addBands(tcDistBG).addBands(tcDistGW).addBands(tcDistBW)

			return img

		img = getTasseledCap(img)
		img = addTCAngles(img)
		return img

	def ND_blue_green(self,img):
		img = img.addBands(img.normalizedDifference(['blue','green']).rename(['ND_blue_green']))
		return img

	def ND_blue_red(self,img):
		img = img.addBands(img.normalizedDifference(['blue','red']).rename(['ND_blue_red']))
		return img

	def ND_blue_nir(self,img):
		img = img.addBands(img.normalizedDifference(['blue','nir']).rename(['ND_blue_nir']))
		return img

	def ND_blue_swir1(self,img):
		img = img.addBands(img.normalizedDifference(['blue','swir1']).rename(['ND_blue_swir1']))
		return img

	def ND_blue_swir2(self,img):
		img = img.addBands(img.normalizedDifference(['blue','swir2']).rename(['ND_blue_swir2']))
		return img

	def ND_green_red(self,img):
		img = img.addBands(img.normalizedDifference(['green','red']).rename(['ND_green_red']))
		return img

	def ND_green_nir(self,img):
		img = img.addBands(img.normalizedDifference(['green','nir']).rename(['ND_green_nir']))  # NDWBI
		return img

	def ND_green_swir1(self,img):
		img = img.addBands(img.normalizedDifference(['green','swir1']).rename(['ND_green_swir1']))  # NDSI, MNDWI
		return img

	def ND_green_swir2(self,img):
		img = img.addBands(img.normalizedDifference(['green','swir2']).rename(['ND_green_swir2']))
		return img

	def ND_red_swir1(self,img):
		img = img.addBands(img.normalizedDifference(['red','swir1']).rename(['ND_red_swir1']))
		return img

	def ND_red_swir2(self,img):
		img = img.addBands(img.normalizedDifference(['red','swir2']).rename(['ND_red_swir2']))
		return img

	def ND_nir_red(self,img):
		img = img.addBands(img.normalizedDifference(['nir','red']).rename(['ND_nir_red']))  # NDVI
		return img

	def ND_nir_swir1(self,img):
		img = img.addBands(img.normalizedDifference(['nir','swir1']).rename(['ND_nir_swir1']))  # NDWI, LSWI, -NDBI
		return img

	def ND_nir_swir2(self,img):
		img = img.addBands(img.normalizedDifference(['nir','swir2']).rename(['ND_nir_swir2'])) # NBR, MNDVI
		return img

	def ND_swir1_swir2(self,img):
		img = img.addBands(img.normalizedDifference(['swir1','swir2']).rename(['ND_swir1_swir2']))
		return img

	def R_swir1_nir(self,img):
		# Add ratios
		img = img.addBands(img.select('swir1').divide(img.select('nir')).rename(['R_swir1_nir']));  # ratio 5/4
		return img

	def R_red_swir1(self,img):
		img = img.addBands(img.select('red').divide(img.select('swir1')).rename(['R_red_swir1']));  # ratio 3/5
		return img

	def EVI(self,img):
		#Add Enhanced Vegetation Index (EVI)
		evi = img.expression(
			'2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
			  'NIR': img.select('nir'),
			  'RED': img.select('red'),
			  'BLUE': img.select('blue')
		  }).float()

		img = img.addBands(evi.rename(['EVI']))

		return img

	def SAVI(self,img):
		# Add Soil Adjust Vegetation Index (SAVI)
		# using L = 0.5;
		savi = img.expression(
			'(NIR - RED) * (1 + 0.5)/(NIR + RED + 0.5)', {
			  'NIR': img.select('nir'),
			  'RED': img.select('red')
		  }).float()
		img = img.addBands(savi.rename(['SAVI']))

		return img

	def IBI(self,img):
		# Add Index-Based Built-Up Index (IBI)
		ibi_a = img.expression(
			'2*SWIR1/(SWIR1 + NIR)', {
			  'SWIR1': img.select('swir1'),
			  'NIR': img.select('nir')
			}).rename(['IBI_A'])


		ibi_b = img.expression(
			'(NIR/(NIR + RED)) + (GREEN/(GREEN + SWIR1))', {
			  'NIR': img.select('nir'),
			  'RED': img.select('red'),
			  'GREEN': img.select('green'),
			  'SWIR1': img.select('swir1')
			}).rename(['IBI_B'])

		ibi_a = ibi_a.addBands(ibi_b)
		ibi = ibi_a.normalizedDifference(['IBI_A','IBI_B'])
		img = img.addBands(ibi.rename(['IBI']))

		return img

	def addTopography(self,img): 
		"""  Function to add 30m SRTM elevation and derived slope, aspect, eastness, and 
		northness to an image. Elevation is in meters, slope is between 0 and 90 deg,
		aspect is between 0 and 359 deg. Eastness and northness are unitless and are
		between -1 and 1. """

		# Import SRTM elevation data
		elevation = ee.Image("USGS/SRTMGL1_003")

		# Calculate slope, aspect, and hillshade
		topo = ee.Algorithms.Terrain(elevation)

		# From aspect (a), calculate eastness (sin a), northness (cos a)
		deg2rad = ee.Number(math.pi).divide(180)
		aspect = topo.select(['aspect'])
		aspect_rad = aspect.multiply(deg2rad)
		eastness = aspect_rad.sin().rename(['eastness']).float()
		northness = aspect_rad.cos().rename(['northness']).float()

		# Add topography bands to image
		topo = topo.select(['elevation','slope','aspect']).addBands(eastness).addBands(northness)
		img = img.addBands(topo)
		return img

	def addJRC(self,img):
		""" Function to add JRC Water layers: 'occurrence', 'change_abs', 
			'change_norm', 'seasonality','transition', 'max_extent' """

		jrcImage = ee.Image("JRC/GSW1_0/GlobalSurfaceWater")

		img = img.addBands(jrcImage.select(['occurrence']).rename(['occurrence']))
		img = img.addBands(jrcImage.select(['change_abs']).rename(['change_abs']))
		img = img.addBands(jrcImage.select(['change_norm']).rename(['change_norm']))
		img = img.addBands(jrcImage.select(['seasonality']).rename(['seasonality']))
		img = img.addBands(jrcImage.select(['transition']).rename(['transition']))
		img = img.addBands(jrcImage.select(['max_extent']).rename(['max_extent']))

		return img


	def getIndices(self,img,covariates):	
		""" add indices to image"""
		# self = indices()
		# no need to add indices that are already there
		# see TODO below, can't use removeDuplicates in .map()
		# indices = self.removeDuplicates(covariates,img.bandNames().getInfo())
		indices = covariates

		for item in indices:
			img = self.functionList[item](img)

		return img

	def removeDuplicates(self,covariateList,bands):
		""" function to remove duplicates, i.e. existing bands do not need to be calculated """
		# TODO: this does not scale to being mappable server side (can't use getInfo in mapped functions)
		# would need to EEify this logic to use with ee.List()'s
		return [elem for elem in covariateList if elem not in bands]

	def renameBands(self,image,prefix):
		""" renames bands with prefix """

		bandnames = image.bandNames()

		def mapBands(band):
			band = ee.String(prefix).cat('_').cat(band)
			return band

		bandnames = bandnames.map(mapBands)

		image = image.rename(bandnames)

		return image

addAllTasselCapIndices(img)

Function to get all tasselCap indices

Source code in src\rlcms\covariates.py 
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
def addAllTasselCapIndices(self,img): 
	""" Function to get all tasselCap indices """

	def getTasseledCap(img):
		"""Function to compute the Tasseled Cap transformation and return an image"""

		coefficients = ee.Array([
			[0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863],
			[-0.2848, -0.2435, -0.5436, 0.7243, 0.0840, -0.1800],
			[0.1509, 0.1973, 0.3279, 0.3406, -0.7112, -0.4572],
			[-0.8242, 0.0849, 0.4392, -0.0580, 0.2012, -0.2768],
			[-0.3280, 0.0549, 0.1075, 0.1855, -0.4357, 0.8085],
			[0.1084, -0.9022, 0.4120, 0.0573, -0.0251, 0.0238]
		])

		bands=ee.List(['blue','green','red','nir','swir1','swir2'])

		# Make an Array Image, with a 1-D Array per pixel.
		arrayImage1D = img.select(bands).toArray()

		# Make an Array Image with a 2-D Array per pixel, 6x1.
		arrayImage2D = arrayImage1D.toArray(1)

		componentsImage = ee.Image(coefficients).matrixMultiply(arrayImage2D).arrayProject([0]).arrayFlatten([['brightness', 'greenness', 'wetness', 'fourth', 'fifth', 'sixth']]).float()

		# Get a multi-band image with TC-named bands.
		return img.addBands(componentsImage);	


	def addTCAngles(img):

		""" Function to add Tasseled Cap angles and distances to an image. Assumes image has bands: 'brightness', 'greenness', and 'wetness'."""

		# Select brightness, greenness, and wetness bands	
		brightness = img.select('brightness')
		greenness = img.select('greenness')
		wetness = img.select('wetness')

		# Calculate Tasseled Cap angles and distances
		tcAngleBG = brightness.atan2(greenness).divide(math.pi).rename(['tcAngleBG'])
		tcAngleGW = greenness.atan2(wetness).divide(math.pi).rename(['tcAngleGW'])
		tcAngleBW = brightness.atan2(wetness).divide(math.pi).rename(['tcAngleBW'])
		tcDistBG = brightness.hypot(greenness).rename(['tcDistBG'])
		tcDistGW = greenness.hypot(wetness).rename(['tcDistGW'])
		tcDistBW = brightness.hypot(wetness).rename(['tcDistBW'])
		img = img.addBands(tcAngleBG).addBands(tcAngleGW).addBands(tcAngleBW).addBands(tcDistBG).addBands(tcDistGW).addBands(tcDistBW)

		return img

	img = getTasseledCap(img)
	img = addTCAngles(img)
	return img

addJRC(img)

Function to add JRC Water layers: 'occurrence', 'change_abs', 'change_norm', 'seasonality','transition', 'max_extent'

Source code in src\rlcms\covariates.py 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
def addJRC(self,img):
	""" Function to add JRC Water layers: 'occurrence', 'change_abs', 
		'change_norm', 'seasonality','transition', 'max_extent' """

	jrcImage = ee.Image("JRC/GSW1_0/GlobalSurfaceWater")

	img = img.addBands(jrcImage.select(['occurrence']).rename(['occurrence']))
	img = img.addBands(jrcImage.select(['change_abs']).rename(['change_abs']))
	img = img.addBands(jrcImage.select(['change_norm']).rename(['change_norm']))
	img = img.addBands(jrcImage.select(['seasonality']).rename(['seasonality']))
	img = img.addBands(jrcImage.select(['transition']).rename(['transition']))
	img = img.addBands(jrcImage.select(['max_extent']).rename(['max_extent']))

	return img

addTopography(img)

Function to add 30m SRTM elevation and derived slope, aspect, eastness, and northness to an image. Elevation is in meters, slope is between 0 and 90 deg, aspect is between 0 and 359 deg. Eastness and northness are unitless and are between -1 and 1.

Source code in src\rlcms\covariates.py 
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
def addTopography(self,img): 
	"""  Function to add 30m SRTM elevation and derived slope, aspect, eastness, and 
	northness to an image. Elevation is in meters, slope is between 0 and 90 deg,
	aspect is between 0 and 359 deg. Eastness and northness are unitless and are
	between -1 and 1. """

	# Import SRTM elevation data
	elevation = ee.Image("USGS/SRTMGL1_003")

	# Calculate slope, aspect, and hillshade
	topo = ee.Algorithms.Terrain(elevation)

	# From aspect (a), calculate eastness (sin a), northness (cos a)
	deg2rad = ee.Number(math.pi).divide(180)
	aspect = topo.select(['aspect'])
	aspect_rad = aspect.multiply(deg2rad)
	eastness = aspect_rad.sin().rename(['eastness']).float()
	northness = aspect_rad.cos().rename(['northness']).float()

	# Add topography bands to image
	topo = topo.select(['elevation','slope','aspect']).addBands(eastness).addBands(northness)
	img = img.addBands(topo)
	return img

getIndices(img, covariates)

add indices to image

Source code in src\rlcms\covariates.py 
242
243
244
245
246
247
248
249
250
251
252
253
def getIndices(self,img,covariates):	
	""" add indices to image"""
	# self = indices()
	# no need to add indices that are already there
	# see TODO below, can't use removeDuplicates in .map()
	# indices = self.removeDuplicates(covariates,img.bandNames().getInfo())
	indices = covariates

	for item in indices:
		img = self.functionList[item](img)

	return img

removeDuplicates(covariateList, bands)

function to remove duplicates, i.e. existing bands do not need to be calculated

Source code in src\rlcms\covariates.py 
255
256
257
258
259
def removeDuplicates(self,covariateList,bands):
	""" function to remove duplicates, i.e. existing bands do not need to be calculated """
	# TODO: this does not scale to being mappable server side (can't use getInfo in mapped functions)
	# would need to EEify this logic to use with ee.List()'s
	return [elem for elem in covariateList if elem not in bands]

renameBands(image, prefix)

renames bands with prefix

Source code in src\rlcms\covariates.py 
261
262
263
264
265
266
267
268
269
270
271
272
273
274
def renameBands(self,image,prefix):
	""" renames bands with prefix """

	bandnames = image.bandNames()

	def mapBands(band):
		band = ee.String(prefix).cat('_').cat(band)
		return band

	bandnames = bandnames.map(mapBands)

	image = image.rename(bandnames)

	return image

returnCovariates(img)

Workflow for computing Landsat and covariates. bands and covariates are hardcoded inside the function.

Source code in src\rlcms\covariates.py 
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
def returnCovariates(img):
	"""Workflow for computing Landsat and covariates. bands and covariates are hardcoded inside the function."""
	# hard coded for now
	bands = ['blue','green','red','nir','swir1', 'swir2']	
	bandLow = ['p20_blue','p20_green','p20_red','p20_nir','p20_swir1', 'p20_swir2']
	bandHigh = ['p80_blue','p80_green','p80_red','p80_nir','p80_swir1', 'p80_swir2']

	"""Calculate the urban, builtup cropland rice and barren primitives """
	covariates = ["ND_blue_green","ND_blue_red","ND_blue_nir","ND_blue_swir1","ND_blue_swir2", \
				  "ND_green_red","ND_green_nir","ND_green_swir1","ND_green_swir2","ND_red_swir1",\
				  "ND_red_swir2","ND_nir_red","ND_nir_swir1","ND_nir_swir2","ND_swir1_swir2",\
				  "R_swir1_nir","R_red_swir1","EVI","SAVI","IBI"]

	index = indices()

	def scaleLandsat(img):
		"""Landast is scaled by factor 0.0001 """
		thermalBand = ee.List(['thermal'])
		thermal = ee.Image(img).select(thermalBand).divide(10)

		otherBands = ee.Image(img).bandNames().removeAll(thermalBand)

		scaled = ee.Image(img).select(otherBands).multiply(0.0001)
		image = ee.Image(scaled.addBands(thermal))        

		return ee.Image(image.copyProperties(img))

	img = scaleLandsat(img)

	def addIndices(img,prefix):
		img = index.addAllTasselCapIndices(img)
		img = index.getIndices(img,covariates)
		img = index.addJRC(img).unmask(0)
		img = index.addTopography(img).unmask(0)	
		if len(prefix) > 0:	
			img = index.renameBands(img,prefix)
		return img


	down = addIndices(img.select(bandLow,bands),"p20")
	middle = addIndices(img.select(bands),"")
	up = addIndices(img.select(bandHigh,bands),"p80")

	img = down.addBands(middle).addBands(up)

	return img

returnCovariatesFromOptions(img, **kwargs)

Computes and adds image covariates according to user settings args: img (ee.Image): image to compute covariates kwargs (dict): a settings dictionary returns: img (ee.Image): multi-band image with all desired covariates

Source code in src\rlcms\covariates.py 
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
def returnCovariatesFromOptions(img,**kwargs):
	"""
	Computes and adds image covariates according to user settings
	args:
		img (ee.Image): image to compute covariates 
		kwargs (dict): a settings dictionary 
	returns:
		img (ee.Image): multi-band image with all desired covariates
	"""
	settings = kwargs
	if 'indices' in kwargs.keys():
		if len(kwargs['indices']) > 0:
			covariates = settings['indices']
			index = indices()

		img = ee.Image(img)
		img = index.getIndices(img,covariates)

	if 'addTasselCap' in kwargs.keys():
		if kwargs['addTasselCap']:
			img = index.addAllTasselCapIndices(img)

	return img