-
Notifications
You must be signed in to change notification settings - Fork 81
/
Copy pathTopographicCCorrection.py
154 lines (119 loc) · 6.04 KB
/
TopographicCCorrection.py
1
2
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
import numpy as np
from numpy import pi
#import datetime
import scipy.stats as stats
#from datetime import timedelta
#import sys
#import os
#import pickle
#debug_logs_directory = r'C:\PROJECTS\gbrunner-raster-functions\pickles'
# Based on QA Band - https://landsat.usgs.gov/collectionqualityband
#QA_BAND_NUM = 7
#misc = [0, 1]
LANDSAT_4_7_CLEAR_PIX_VALS = [672, 676, 680, 684]
LANDSAT_8_CLEAR_PIX_VALS = [20480, 20484, 20512, 23552]#[2720, 2724, 2728, 2732]
LANDSAT_CLEAR_PIX_VALS = LANDSAT_4_7_CLEAR_PIX_VALS + LANDSAT_8_CLEAR_PIX_VALS
# Terrain C-Correction based on the following papers:
# 1. Teillet, Guindon, and Goodenough (1982)
# 2. Sola et al (2016)
# 3. Holben et al (1980) https://www.asprs.org/wp-content/uploads/pers/1980journal/sep/1980_sep_1191-1200.pdf
# 4. Dr. Mort Canty - https://github.com/mortcanty/CRCPython
class TopographicCCorrection():
def __init__(self):
self.name = 'Topographic C Correction'
self.description = 'Topographic C-Correction based on the paper from Teillet, Guindon, and Goodenough (1982).'
self.metadata = []
def getParameterInfo(self):
return [
{
'name': 'rasters',
'dataType': 'rasters',
'value': None,
'required': True,
'displayName': 'Rasters',
'description': 'The collection of overlapping rasters to aggregate.'
},
{
'name': 'slope',
'dataType': 'raster',
'value': None,
'required': True,
'displayName': "Slope",
'description': "Slope Derived from Digital Elevation Model."
},
{
'name': 'aspect',
'dataType': 'raster',
'value': None,
'required': True,
'displayName': "Aspect",
'description': "Aspect Derived from Digital Elevation Model."
}
]
def getConfiguration(self, **scalars):
return {
'inheritProperties': 4 | 8, # inherit everything but the pixel type (1) and NoData (2)
'invalidateProperties': 2 | 4, # invalidate histogram and statistics because we are modifying pixel values
'inputMask': True, # need raster mask of all input rasters in .updatePixels().
'resampling': False, # process at native resolution
'keyMetadata': ['SunAzimuth','SunElevation','AcquisitionDate']
}
def updateRasterInfo(self, **kwargs):
#outStats = {'minimum': -1, 'maximum': 1}
#self.outBandCount = 6
kwargs['output_info']['pixelType'] = 'f4' # output pixels are floating-point values
kwargs['output_info']['histogram'] = () # no statistics/histogram for output raster specified
kwargs['output_info']['statistics'] = () # outStatsTuple
self.metadata = kwargs['rasters_keyMetadata']
return kwargs
def updateKeyMetadata(self, names, bandIndex, **keyMetadata):
return keyMetadata
def updatePixels(self, tlc, shape, props, **pixelBlocks):
#fname = '{:%Y_%b_%d_%H_%M_%S}_t.txt'.format(datetime.datetime.now())
#filename = os.path.join(debug_logs_directory, fname)
image_pix_blocks = pixelBlocks['rasters_pixels']
image_pix_array = np.asarray(image_pix_blocks)
##pickle_filename = os.path.join(debug_logs_directory, fname)
##pickle.dump(image_pix_array, open(pickle_filename[:-4]+'landsat.p',"wb"))
slope_pix_blocks = pixelBlocks['slope_pixels']
slope_pix_array = np.asarray(slope_pix_blocks)
slope_rads = slope_pix_array * pi/180
##pickle_filename = os.path.join(debug_logs_directory, fname)
##pickle.dump(slope_rads, open(pickle_filename[:-4]+'slope.p',"wb"))
aspect_pix_blocks = pixelBlocks['aspect_pixels']
aspect_pix_array = np.asarray(aspect_pix_blocks)
aspect_rads = aspect_pix_array * pi/180
##pickle_filename = os.path.join(debug_logs_directory, fname)
##pickle.dump(aspect_rads, open(pickle_filename[:-4]+'aspect.p',"wb"))
t_vals = [j['acquisitiondate'] for j in self.metadata]
sun_az = [j['sunazimuth'] for j in self.metadata]
sun_el = [j['sunelevation'] for j in self.metadata]
#https://en.wikipedia.org/wiki/Solar_zenith_angle
sun_ze = [(90 - el) for el in sun_el]
#Equation
#Corrected Image = Image * (cos(solar zenith angle) + C)/(cos(solar incidence angle) + C)
#Where C is an empiraicle parameter
#cos(solar incidence angle) =
pix_array_dim = image_pix_array.shape
num_bands = pix_array_dim[1]
num_squares_x = pix_array_dim[2]
num_squares_y = pix_array_dim[3]
result = np.zeros((num_bands, num_squares_x, num_squares_y))
if len(sun_az) == 1:
sun_az_rad = sun_az[0] * pi / 180
sun_ze_rad = sun_ze[0] * pi / 180
# Topographic Effect on Spectral Response from Nadir-Pointing Sensors
# https://www.asprs.org/wp-content/uploads/pers/1980journal/sep/1980_sep_1191-1200.pdf
cos_i = (np.cos(slope_rads[0, :, :]) * np.cos(sun_ze_rad)) + \
(np.sin(slope_rads[0, :, :] * np.sin(sun_ze_rad)) * (np.cos(sun_az_rad - aspect_rads[0, :, :])))
for k in range(num_bands):
image_flattened = image_pix_array[0, k, :, :].ravel()
cos_gamma = cos_i.ravel()
# Adding linear regression from Dr. Mort Canty - https://github.com/mortcanty/CRCPython
m, b, r, _, _ = stats.linregress(cos_gamma, image_flattened)
C = b / m
result[k, :, :] = image_pix_array[0, k, :, :] * (np.cos(sun_ze_rad) + C) / (cos_i + C)
mask = np.ones((pix_array_dim[1], num_squares_x, num_squares_y))
pixelBlocks['output_mask'] = mask.astype('u1', copy = False)
pixelBlocks['output_pixels'] = result.astype(props['pixelType'], copy=False)
return pixelBlocks