-
Notifications
You must be signed in to change notification settings - Fork 15
/
MODIS_TFA_Daily_Mean_Export_2_Assets
242 lines (212 loc) · 9.81 KB
/
MODIS_TFA_Daily_Mean_Export_2_Assets
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
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
// This code uses Temporal Fourier Analysis to build 365 images, for each day of year,
// of the climatological tempratures based on the daily mean temperatures.
// Change the path to your own assets folder
var Assets_path = 'users/shilosh/MODIS_TFA_Daily_'
var Temperature_Band = 'LST_Daily_Mean';
var Day_Temperature_Band = 'LST_Day_1km';
var Night_Temperature_Band = 'LST_Night_1km';
var collection = 'MODIS/006/MYD11A1';
// var geometry = ee.Geometry.Rectangle([33.2,29.0,36.6,34.0]); //Israel
var geometry = ee.Geometry.Rectangle([20,50,25,54]); //Poland
// var geometry = ee.Geometry.Rectangle([43,0,48,4]); //France
// var geometry = ee.Geometry.Rectangle([49,7,54,11]); //Germany
// var geometry = ee.Geometry.Rectangle([-76,36.0,-79,40]); //Maryland
// var geometry = ee.Geometry.Rectangle([-90,41,-95,44.0]);//Iowa
// var geometry = ee.Geometry.Rectangle([-60,-6,-55,-2]); //Brazil
// var geometry = ee.Geometry.Rectangle([-50,-20,-45,-16]); //SouthBrazil
// var geometry = ee.Geometry.Rectangle([-69,-34,-64,-30]); //Argentina
// var geometry = ee.Geometry.Rectangle([12,24,17,28]); //Libya (Sahara)
// var geometry = ee.Geometry.Rectangle([105,30,110,34]); //China
// var geometry = ee.Geometry.Rectangle([145,-37,150,-33]); //SouthAustralia
// var geometry = ee.Geometry.Rectangle([135,-25,140,-20]); //MidAustralia
// var geometry = ee.Geometry.Rectangle([-130,58,-125,62]); //Canada (Mid latitudes)
// var TFA_File_Name = 'MODIS_Daily_TFA_Israel'
var TFA_File_Name = 'MODIS_Daily_TFA_Poland'
// var TFA_File_Name = 'MODIS_Daily_TFA_France'
// var TFA_File_Name = 'MODIS_Daily_TFA_Germany'
// var TFA_File_Name = 'MODIS_Daily_TFA_Maryland'
// var TFA_File_Name = 'MODIS_Daily_TFA_Iowa'
// var TFA_File_Name = 'MODIS_Daily_TFA_Brazil'
// var TFA_File_Name = 'MODIS_Daily_TFA_SouthBrazil'
// var TFA_File_Name = 'MODIS_Daily_TFA_Argentina'
// var TFA_File_Name = 'MODIS_Daily_TFA_Libya'
// var TFA_File_Name = 'MODIS_Daily_TFA_China'
// var TFA_File_Name = 'MODIS_Daily_TFA_SouthAustralia'
// var TFA_File_Name = 'MODIS_Daily_TFA_MidAustralia'
// var TFA_File_Name = 'MODIS_Daily_TFA_Canada'
var geometry_json = geometry.toGeoJSON();
// The number of cycles per year to model.
var harmonics = 3;
// Make a list of harmonic frequencies to model.
// These also serve as band name suffixes.
var harmonicFrequencies = ee.List.sequence(1, harmonics);
// Function to get a sequence of band names for harmonic terms.
var getNames = function(base, list) {
return ee.List(list).map(function(i) {
return ee.String(base).cat(ee.Number(i).int());
});
};
// Construct lists of names for the harmonic terms.
var cosNames = getNames('cos_', harmonicFrequencies);
var sinNames = getNames('sin_', harmonicFrequencies);
//convert Kelvin to Celsius
var k2celsius = function(image) {
return image.updateMask(image.select(Day_Temperature_Band))
.updateMask(image.select(Night_Temperature_Band))
.reduce( ee.Reducer.mean()).rename(Temperature_Band)
.multiply(ee.Image(0.02))
.subtract(ee.Image(273.15)) //convert Kelvin to Celsius
.set('doy', image.get('doy'));
};
// Add a band with doy to the colection.
function createDoyBand(img) {
var d = ee.Number(img.get('doy'))
return img.addBands(ee.Image(d).rename('t'));
}
// Function to compute the specified number of harmonics
// and add them as bands.
// for p = 0, Nharmonics-1 do begin
// a[p] = a[p] + 2*TS_t*cos(2*!pi*(p+1)*(t+1)/N)
// b[p] = b[p] + 2*TS_t*sin(2*!pi*(p+1)*(t+1)/N)
var addHarmonics = function(freqs) {
return function(image) {
// Make an image of frequencies.
var frequencies = ee.Image.constant(freqs);
// This band should represent lst.
var lst = ee.Image(image).select(Temperature_Band);
// This band should represent days from start.
var time = ee.Image(image).select('t');
// Get the cosine terms.
var cosines = lst.multiply(2.0).multiply((frequencies.multiply(2.0 * Math.PI).multiply(time).divide(365.0)).cos())
.rename(cosNames);
// Get the sin terms.
var sines = lst.multiply(2.0).multiply((frequencies.multiply(2.0 * Math.PI).multiply(time).divide(365.0)).sin())
.rename(sinNames);
return image.addBands(cosines).addBands(sines);
};
};
var LST_Day = ee.ImageCollection(ee.List.sequence(1, 365).map(function (doy){
return ee.ImageCollection(collection)
.select(Day_Temperature_Band)
.filter(ee.Filter.calendarRange(doy, doy, 'day_of_year'))
.mean()
.set({doy: doy})
//.addBands(ee.Image(ee.Number(doy)).rename('t'))
}))
var LST_Night = ee.ImageCollection(ee.List.sequence(1, 365).map(function (doy){
return ee.ImageCollection(collection)
.select(Night_Temperature_Band)
.filter(ee.Filter.calendarRange(doy, doy, 'day_of_year'))
.mean()
.set({doy: doy})
//.addBands(ee.Image(ee.Number(doy)).rename('t'))
}))
print('LST_Day = ', LST_Day)
// // For each day cloudy pixels with no data, insert yearly mean values.
// var LST_Mean
// function BlendWithYearlyMean(img) {
// return img.blend(LST_Mean);
// }
// LST_Mean = LST_Day.mean()
// LST_Day = LST_Day.map(BlendWithYearlyMean)
// LST_Mean = LST_Night.mean()
// LST_Night = LST_Night.map(BlendWithYearlyMean)
// Use an equals filter to specify how the collections match.
var Filter = ee.Filter.equals({
leftField: 'doy',
rightField: 'doy'
});
// Define the join.
var innerJoin = ee.Join.inner('primary', 'secondary');
// Join LST_Day with LST_Night by DOY
// Apply the join.
var CFSV2_JoinInner = innerJoin.apply(LST_Day, LST_Night, Filter);
// Add band
var LST = CFSV2_JoinInner.map(function(f) {
var Night = ee.Image(f.get('secondary'));
var Day = ee.Image(f.get('primary'));
return Day.addBands(Night);
})
LST = ee.ImageCollection(LST)
.map(k2celsius)
.map(createDoyBand)
print('LST = ', LST)
LST = LST.map(addHarmonics(harmonicFrequencies));
// a = a / N_TS
// b = b / N_TS
var num_days = 365;//ee.Date(lastDay).difference(ee.Date(firstDay), 'day');
var a_coef_1 = LST.select('cos_1').sum().divide(num_days);
var a_coef_2 = LST.select('cos_2').sum().divide(num_days);
var a_coef_3 = LST.select('cos_3').sum().divide(num_days);
var b_coef_1 = LST.select('sin_1').sum().divide(num_days);
var b_coef_2 = LST.select('sin_2').sum().divide(num_days);
var b_coef_3 = LST.select('sin_3').sum().divide(num_days);
// ; ========================================
// ; prepare the first 3 harmonics of the TFA
// ; ========================================
// H0 = total(TS) / N_TS ;mean
// FTA = MAKE_ARRAY(N_TS, /FLOAT, VALUE=0.)
// H = MAKE_ARRAY(Nharmonics, N_TS, /FLOAT, VALUE=0.)
// omega = MAKE_ARRAY(Nharmonics, /FLOAT, VALUE=0.)
// Phase = MAKE_ARRAY(Nharmonics, /FLOAT, VALUE=0.)
// Amplitude = MAKE_ARRAY(Nharmonics, /FLOAT, VALUE=0.)
// for p = 0, Nharmonics-1 do begin
// omega[p] = 2. * !pi * (p+1) / N
// Phase[p] = atan( -b[p], a[p] )
// Amplitude[p] = sqrt( a[p]^2 + b[p]^2 )
var H0 = LST.select(Temperature_Band).mean();
// Get the omega terms.
var omega_1 = ee.Image(1.0).multiply(2.0 * Math.PI).divide(365.0);
var omega_2 = ee.Image(2.0).multiply(2.0 * Math.PI).divide(365.0);
var omega_3 = ee.Image(3.0).multiply(2.0 * Math.PI).divide(365.0);
var phase_1 = a_coef_1.atan2(b_coef_1.multiply(ee.Image(-1.0)));
var phase_2 = a_coef_2.atan2(b_coef_2.multiply(ee.Image(-1.0)));
var phase_3 = a_coef_3.atan2(b_coef_3.multiply(ee.Image(-1.0)));
var amplitude_1 = (a_coef_1.pow(2).add(b_coef_1.pow(2))).sqrt();
var amplitude_2 = (a_coef_2.pow(2).add(b_coef_2.pow(2))).sqrt();
var amplitude_3 = (a_coef_3.pow(2).add(b_coef_3.pow(2))).sqrt();
// for t = 0, N_TS-1 do begin
// H[p,t] = Amplitude[p] * cos(omega[p]*(t+1) + Phase[p])
// Function to add a H band.
var addH = function(image) {
var H_1 = amplitude_1.multiply((omega_1.multiply(image.select('t')).add(phase_1)).cos()).rename('H_1');
var H_2 = amplitude_2.multiply((omega_2.multiply(image.select('t')).add(phase_2)).cos()).rename('H_2');
var H_3 = amplitude_3.multiply((omega_3.multiply(image.select('t')).add(phase_3)).cos()).rename('H_3');
// FTA[t] = H0 + H[0,t] + H[1,t] + H[2,t]
var TFA = H0.add(H_1).add(H_2).add(H_3).rename('TFA');
var Anomalies = TFA.subtract(image.select(Temperature_Band)).rename('Anomalies');
return image.addBands(TFA).addBands(Anomalies);
};
var LST_TFA = LST.map(addH);
// Define the visualization parameters.
var vizParams = {
bands: ['mean', 'phase', 'amplitude'],
min: [285.0,2.0,5.0],
max: [295.0,3.0,10.0],
gamma: [1.0, 1.0, 1.0]
};
//display the image.
var img = ee.Image(H0.addBands(phase_1).addBands(amplitude_1).rename('mean', 'phase', 'amplitude'));
// Map.addLayer(img, vizParams, 'RGB');
// var modelTFA = LST_TFA.select(['TFA', Temperature_Band, 'Anomalies']);
// var modisTFA = LST_TFA.select(['TFA', Temperature_Band]);
// print(LST_TFA)
//Iterating over the image collection using this function....
// var geometry = ee.FeatureCollection('ft:19Un8ogYFYTnR9p_VsMQ1tlzbR-bVwZoj7ipKy-bC');
var TFA_Images = LST_TFA.select('TFA')
.sort('t')
.iterate(function(img, all) {
return ee.Image(all).addBands(img);
}, ee.Image().select());
// print('TFA_Images = ', TFA_Images );
// Map.addLayer(ee.Image(TFA_Images))
// Export one single file with all the bands
Export.image.toAsset({
image: ee.Image(TFA_Images),
assetId: Assets_path + TFA_File_Name,
region: geometry_json,
crs: ee.Image('MODIS/006/MYD11A1/2018_01_01').projection().crs().getInfo(),
crsTransform: [926.625433056,0,-20015109.354,0,-926.625433055,10007554.677],
description: 'exportToAsset-MODIS-TFA-' + TFA_File_Name,
pyramidingPolicy: {'.default': 'sample'}
});