Vegetation Indices#
MODIS NDVI#
Data download from Google Earth Engine#
The MOD13Q1 V6 product provides a Vegetation Index (VI) value at a per pixel basis. There are two primary vegetation layers. The first is the Normalized Difference Vegetation Index (NDVI) which is referred to as the continuity index to the existing National Oceanic and Atmospheric Administration-Advanced Very High Resolution Radiometer (NOAA-AVHRR) derived NDVI. The second vegetation layer is the Enhanced Vegetation Index (EVI) that minimizes canopy background variations and maintains sensitivity over dense vegetation conditions. The MODIS NDVI and EVI products are computed from atmospherically corrected bi-directional surface reflectances that have been masked for water, clouds, heavy aerosols, and cloud shadows.
To compute and download the mean anual MODIS NDVI data from google earth engine. Open the Google Earth Engine Code and paste the lines of code provided below
1 /**** Start of imports. If edited, may not auto-convert in the playground. ****/
2 var table = ee.FeatureCollection("users/derickongeri/NorthAfrica");
3 /***** End of imports. If edited, may not auto-convert in the playground. *****/
4 var dataset = ee.ImageCollection('MODIS/006/MOD13Q1')
5 .filter(ee.Filter.date('2018-01-01', '2019-01-01'))
6 .mean()
7 .clip(table);
8
9 var ndvi = dataset.select('NDVI');
10 var ndviVis = {
11 min: 0.0,
12 max: 8000.0,
13 palette: [
14 'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
15 '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
16 '012E01', '011D01', '011301'
17 ],
18 };
19 Map.centerObject(table);
20 Map.addLayer(ndvi, ndviVis, 'NDVI');
21
22 Export.image.toCloudStorage({
23 image:ndvi,
24 description: 'NDVI',
25 maxPixels:1e13,
26 scale:250,
27 bucket:'oss_ldms_vi',
28 region:table
29 })
Landsat derived NDVI, MSAVI, SAVI#
Data download from Google Earth Engine#
Several annual NDVI composite techniques have been discussed to overcome current Landsat 5 artifacts. MISLAND uses annual NDVI products based on different percentiles in order to better qualify and quantify the Vegetation Loss Index.
To compute and download desired percentile composites from google earth engine. Open the Google Earth Engine Code and paste the lines of code provided below
1 // Required Data Inputs
2 // ===================
3 // * USGS/NASA's Landsat 4 surface reflectance tier 1 dataset (August 1982 - December 1993)
4 // * USGS/NASA's Landsat 5 surface reflectance tier 1 dataset (January 1, 1984 - May 5, 2012)
5 // * USGS/NASA's Landsat 7 surface reflectance tier 1 dataset (January 1, 1999 - December 31, 2019)
6 // * USGS/NASA's Landsat 8 surface reflectance tier 1 dataset (April 11, 2013 - December 31, 2019)
7 // * Study Area Polygon
8
9 var countries = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017");
10
11 var Year='2002'
12
13
14 var country = 'Libya';
15 var ALGO = "PERCENTILE65"; // PERCENTILE75 // PERCENTILE65 // PERCENTILE60 // MEDIAN
16 var cloudCoveragePercentage = 80;
17
18
19 var studyArea = countries.filter(ee.Filter.eq('country_na',country ))
20 //Map.addLayer(studyArea);
21
22 /*
23 borders are quite coarse
24 var northAfrica = ee.FeatureCollection('users/derickongeri/Admin')
25 // Replace country name with EGYPT, LIBYA, ALGERIA, MAURITANIA, MOROCCO, TUNISIA
26 var country = 'TUNISIA';
27 var studyArea = northAfrica.filter(ee.Filter.eq('NAME', country))
28 Map.addLayer(studyArea);
29 */
30
31
32
33
34 var start_date = Year+ '-01-01';
35 var end_date = Year+ '-12-31';
36
37
38 //--------------------------------------------------------------------
39 // Landsat 4, 5, 7 cloudmask
40 //--------------------------------------------------------------------
41
42 // If the cloud bit (5) is set and the cloud confidence (7) is high
43 // or the cloud shadow bit is set (3), then it's a bad pixel.
44 var cloudMaskL7 = function(image) {
45 var qa = image.select('pixel_qa');
46 var cloud = qa.bitwiseAnd(1 << 5)
47 .and(qa.bitwiseAnd(1 << 7))
48 .or(qa.bitwiseAnd(1 << 3));
49
50 // Remove edge pixels that don't occur in all bands
51 //var mask2 = image.mask().reduce(ee.Reducer.min())//.focal_min(300,'square','meters').eq(0);
52 //var mask2 = image.select('B4').reduce(ee.Reducer.min()).gt(0)//.focal_min(500,'square','meters');
53 // Remove edge pixels that don't occur in all bands
54 var mask3 =
55 (image.select('B3').gt(100))
56 .and(image.select('B4').gt(100))
57
58
59 .and(image.select('B4').lt(10000))
60 .and(image.select('B3').lt(10000))
61
62
63
64 return image.updateMask(cloud.not()).updateMask(mask3)//.updateMask(mask2)//.clip(image.geometry().buffer(-5000))//.or(mask3));
65 };
66
67
68 var cloudMaskL45 = function(image) {
69 var qa = image.select('pixel_qa');
70 var cloud = qa.bitwiseAnd(1 << 5)
71 .and(qa.bitwiseAnd(1 << 7))
72 .or(qa.bitwiseAnd(1 << 3));
73
74 // Remove edge pixels that don't occur in all bands
75 //var mask2 = image.mask().reduce(ee.Reducer.min());
76 var mask2 =
77 (image.select('B3').gt(100))
78 .and(image.select('B4').gt(100))
79
80
81 .and(image.select('B4').lt(10000))
82 .and(image.select('B3').lt(10000))
83
84 return (image.updateMask(cloud.not()).updateMask(mask2))//.clip(image.geometry().buffer(-5000))//.updateMask(mask2);
85 };
86
87
88 //--------------------------------------------------------------------
89 // Landsat 8 cloudmask
90 //--------------------------------------------------------------------
91
92 // Bits 3 and 5 are cloud shadow and cloud, respectively.
93 function maskL8sr(image) {
94 var cloudShadowBitMask = (1 << 3);
95 var cloudsBitMask = (1 << 5);
96
97 // Get the pixel QA band.
98 var qa = image.select('pixel_qa');
99
100 // Both flags should be set to zero, indicating clear conditions.
101 var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
102 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
103 var mask2 =
104
105 (image.select('B5').gt(100))
106 .and(image.select('B4').gt(100))
107
108
109 .and(image.select('B5').lt(10000))
110 .and(image.select('B4').lt(10000))
111
112 //var mask2 = image.mask().reduce(ee.Reducer.min()).focal_min(500,'square','meters');
113 //return image
114 return image.updateMask(mask).updateMask(mask2)//.clip(image.geometry().buffer(-5000));
115 }
116
117
118
119
120 // Apply Cloudmask to L4.5.7
121 var L4 = ee.ImageCollection("LANDSAT/LT04/C01/T1_SR")
122 .filterDate(start_date, end_date)
123 .filter(ee.Filter.lessThan('CLOUD_COVER_LAND', cloudCoveragePercentage))
124 .filterBounds(studyArea)
125 .map(cloudMaskL45)
126 .select(['B3', 'B4'], ['RED', 'NIR']);;
127
128 var L5 = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
129 .filterDate(start_date, end_date)
130 .filter(ee.Filter.lessThan('CLOUD_COVER_LAND', cloudCoveragePercentage))
131 .filterBounds(studyArea)
132 .map(cloudMaskL45)
133 .select(['B3', 'B4'], ['RED', 'NIR']);;
134
135 var L7a = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
136 .filterDate('1999-01-01', '2003-04-01')
137 .filterDate(start_date, end_date)
138 .filter(ee.Filter.lessThan('CLOUD_COVER_LAND', 100))
139 .filterBounds(studyArea)
140 .map(cloudMaskL7)
141 .select(['B3', 'B4'], ['RED', 'NIR']);;
142 var L7b = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
143 .filterDate('2012-01-01', '2013-12-31')
144 .filterDate(start_date, end_date)
145 .filter(ee.Filter.lessThan('CLOUD_COVER_LAND', 100))
146 .filterBounds(studyArea)
147 .map(cloudMaskL7)
148 .select(['B3', 'B4'], ['RED', 'NIR']);;
149
150 var L7 = L7a.merge(L7b);
151
152 var L8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
153 .filterDate(start_date, end_date)
154 .filter(ee.Filter.lessThan('CLOUD_COVER', cloudCoveragePercentage))
155 .filterBounds(studyArea)
156 //.filterBounds(AOI)
157 .map(maskL8sr)
158 .select(['B4', 'B5'], ['RED', 'NIR']);
159
160
161
162
163
164
165 //Define collection
166
167
168 //--------------------------------------------------------------------
169 // Merge Landsat 4, 5, 8 imagery collections and filter all by date/place
170 //--------------------------------------------------------------------
171
172 //Merge Landsat 4, 5 , 7 '
173
174 var L4578 = L4.merge(L5).merge(L7).merge(L8);
175
176
177 //--------------------------------------------------------------------
178 // Create NDVI Collection
179 //--------------------------------------------------------------------
180
181 var NDVI = function(image) {
182 return image.normalizedDifference(['NIR', 'RED']).rename('NDVI');
183 //return image.addBands(ndvi);
184 };
185
186 if (ALGO=='MEDIAN'){
187 var suffix = 'median';
188 var annualNDVI = L4578.map(NDVI).median().clip(studyArea);
189 }
190
191 if (ALGO=='PERCENTILE75'){
192 var suffix = '75pc';
193 var annualNDVI = L4578.map(NDVI).reduce(ee.Reducer.percentile([75])).clip(studyArea);
194 }
195
196 if (ALGO=='PERCENTILE65'){
197 var suffix = '65pc';
198 var annualNDVI = L4578.map(NDVI).reduce(ee.Reducer.percentile([65])).clip(studyArea);
199 }
200
201
202 var ndvi_visualization = {
203 min: -0.22789797020331423,
204 max: 0.6575894075894075,
205 palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
206 '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
207 };
208 Map.addLayer(annualNDVI, ndvi_visualization, 'NDVI');
209
210 //--------------------------------------------------------------------
211 // Export as GeoTIFF
212 //--------------------------------------------------------------------
213
214
215 Export.image.toDrive({
216 image: annualNDVI,
217 description: country + '_NDVI_' + suffix + '_' + Year,
218 scale: 30,
219 region: studyArea,
220 maxPixels: 1e13,
221 fileFormat: 'GeoTIFF',
222 folder:'GEE_classification',
223 formatOptions: {
224 cloudOptimized: true
225 },
226 skipEmptyTiles: true
227 });
228
229 //Map.addLayer(L7.first(), {}, 'L7');