Sahara-Cloudiness-with-MODIS-and-GEE_js

summer 2024 by Miha

Study area

~ 5 000 000 km2 (little less)

Clouds area retrieving method:

Imagery analytics where performed on Google Earth Engine(GEE) platform using java script and MODIS Surface Reflectance product: MOD09GA version 6.1;
NDCI = Normalized Difference Cloud Index

NDCI = (sur_refl_b01 - sur_refl_b06)/(sur_refl_b01 + sur_refl_b06).

Imagery processing:
1 SDS(C_flag)State_1km Bitmask QA; cloud state bit was used as a mask for cloud mapping;
2 LOW PASS includes pixels with NDCI greather than -0.14
3 MID PASS processes NDCI greather than 0
4 HIGH PASS processes red band(b01) surface reflectance greather than 0.6.

_{*fast grafic presentation: screenshoots from GEE (~)center is at three borders(Algeria,Mali,Niger); RGB(MODISTerra SR 10.jan.2015);}

WHITE = QA(SDS cloud flag)
BLUE = low_pass
YELLOW = mid_pass
RED = high_pass

_{Used MODIS bands:
500m Surface Reflectance Band 1: (620-670 nm); red
500m Surface Reflectance Band 6: (1628-1652 nm); SWIR I
1000m state_1km: Reflectance data state QA flags, SDS(Scientific Data Set), cloud state}

Results:

Average cloudiness of study area is ~ 8 %. The cloudiest year was 2016, 9.8% and year 2001 was the least cloudy year with average cloudiness 5.6%.

Main Java script for retrieving data

Script is optimized to process and export results from large data input where output preserves input bands resolution.There is no GUI, just table export.

/*
 * @license
 * Author: Miha Žemva
 * correspondence: [email protected]
 * Licensed under the Apache License, Version 2.0 
 * 
 *     https://www.apache.org/licenses/LICENSE-2.0
*/
var MODIS_SR_coll = ee.ImageCollection('MODIS/061/MOD09GA') //MODIS/061/MYD09GA
       .filterDate('2000-03-01', '2000-04-01')
       .select(['state_1km','sur_refl_b01','sur_refl_b06'])
       .map(function(i){return i.clip(Sahara_study)});
       
var low_pass_th = -0.14;
var mid_pass_th = 0;
var red_high_pass = 6000;

function Clouds_MOD09GA(i) {
    var date = i.date().format('YYYY-MM-dd');
    var year = i.date().get('year');
    var month = i.date().get('month');
    var week_Y = i.date().get('week');
    var day_Y = i.date().format('DDD');
    var C_flag = i.select('state_1km').bitwiseAnd(1 << 0).eq(1).rename('Cflag_state_1km');
    
    var NDCI = i.normalizedDifference(['sur_refl_b01', 'sur_refl_b06']).rename('NDCI_redswir1');
    var low_pass_NDCI = NDCI.gt(low_pass_th).and(C_flag.eq(1)).rename('NDCI_low_pass');
    var mid_pass_NDCI = NDCI.gt(mid_pass_th).and(C_flag.eq(1)).rename('NDCI_mid_pass');
    var high_pass_RED = mid_pass_NDCI.eq(1).and(i.select('sur_refl_b01').gt(red_high_pass)).rename('RED_high_pass');

  var low_pass_C = low_pass_NDCI.reduceRegion({
  reducer: ee.Reducer.mean().unweighted(),
  geometry: Sahara_study,
  scale: 500, 
  maxPixels: 1e9,}).values();
  
  var mid_pass_C = mid_pass_NDCI.reduceRegion({
  reducer: ee.Reducer.mean().unweighted(),
  geometry: Sahara_study,
  scale: 500, 
  maxPixels: 1e9,}).values();
  
  var high_pass_C = high_pass_RED.reduceRegion({
  reducer: ee.Reducer.mean().unweighted(),
  geometry: Sahara_study,
  scale: 500, 
  maxPixels: 1e9,}).values();
  
  var count_GA = i.select('sur_refl_b01').reduceRegion({
  reducer: ee.Reducer.count().unweighted(),
  geometry: Sahara_study,
  scale: 500, 
  maxPixels: 1e9,}).values();
  
  var row = ee.List([date, year, month, week_Y, day_Y, low_pass_C, mid_pass_C, high_pass_C, count_GA]).flatten();

return ee.Feature(null, {'rows1': row});
}

Export.table.toDrive({
  collection: MODIS_SR_coll.map(Clouds_MOD09GA),
  description: 'description',
  folder: 'GEE_sahra',
  fileFormat: 'CSV',
  selectors: ['rows1']
});

MAX cloudiness distribution

MIN cloudiness distribution

_{Study area had zero cloud cover on three days: 20.dec 2006, 7.feb 2010 and 11.mar 2018.}

filtering based on counting valid pixels in study area

Analysed imagery has time interval 1.march 2000-27.feb 2023, that is total 8400 days(68 days had no data on study area; additional 138 days where removed later in Google Sheets, where days with count_GA higher than 18.545.000 px are considered valid.

No. of daily imagery used in analytics

MODIS Terra maps:

Number of cloudy days in one year in each pixel (20 years averages 2000-2020). This maps where exported from GEE and processed in SAGA-GIS; they don't exclude days with significant sensor failure(like px_count in study area), therefore the true values [days/year] are slightly higher(export from GEE was sum). This maps were made before the study area was drawn.

High pass _{QA-C_flag & red_band refl gt(0.6)}

MID pass _{QA-C_flag & NDCI gt(0)}

HIGH pass | MID pass _{Avg no. of days on areas... Areas where generated by reclasifiying, resampling and segmentation of QA_SDS(Cloud flag) map in SAGA-GIS.}

LOW pass _{QA_C_flag & NDCI gt(-0.14)}

QA-C_flag _{state_1km Bitmask (bit 0 = 1)}

LOW pass | QA(C_flag) _{Avg no. of days on areas...}

Visualization Java script

/*
 * @license
 * Author: Miha Žemva
 * correspondence: [email protected]
 * Licensed under the Apache License, Version 2.0 
 * 
 *     https://www.apache.org/licenses/LICENSE-2.0
*/

var MODIS_SR_coll = ee.ImageCollection('MODIS/061/MOD09GA')
       .filterDate('2000-03-01', '2000-03-16')
       .select(['state_1km','sur_refl_b01','sur_refl_b06'])
       .map(function(i){return i.clip(Sahara_study)});

var low_pass_th = -0.14;
var mid_pass_th = 0;
var red_high_pass = 6000;

function Clouds_MOD09GA(i) {
    var C_flag = i.select('state_1km').bitwiseAnd(1 << 0).eq(1).rename('Cflag_state_1km');
    var NDCI = i.normalizedDifference(['sur_refl_b01', 'sur_refl_b06']).rename('NDCI_redswir1');
    var low_pass_NDCI = NDCI.gt(low_pass_th).and(C_flag.eq(1)).rename('NDCI_low_pass');
    var mid_pass_NDCI = NDCI.gt(mid_pass_th).and(C_flag.eq(1)).rename('NDCI_mid_pass');
    var high_pass_RED = mid_pass_NDCI.eq(1).and(i.select('sur_refl_b01').gt(red_high_pass)).rename('RED_high_pass');
  
  return ee.Image([C_flag, low_pass_NDCI, mid_pass_NDCI, high_pass_RED]).copyProperties(i,['system:time_start']);
}

var OUTPUT_images = MODIS_SR_coll.map(Clouds_MOD09GA);

var OUTPUT_List = ee.Image(OUTPUT_images.toList(OUTPUT_images.size()));
print('OUTPUT_List', OUTPUT_List);

var listed_Day = 0;

var OUTPUT_Day = ee.Image(OUTPUT_images.toList(OUTPUT_images.size()).get(listed_Day));
print('listed Day :'+listed_Day, ee.Date((OUTPUT_Day.get('system:time_start'))));
print('Study_area_m2', Sahara_study.area());

var Cmask_chart = ui.Chart.image.seriesByRegion(OUTPUT_images, Sahara_study, 
      ee.Reducer.mean(), 'Cflag_state_1km', 1000, 'system:time_start','label')
        .setOptions({
          title: 'Cloudines MODIS_sr Cmask',
          vAxis: {title: 'band: state_1km Bitmask'},
          pointSize: 2,
          colors: ['E37D05']
        });
print(Cmask_chart);

var NDCI_low_chart = ui.Chart.image.seriesByRegion(OUTPUT_images, Sahara_study, 
      ee.Reducer.mean(),'NDCI_low_pass', 500, 'system:time_start','label')
        .setOptions({
          title: 'Cloudines MODIS_sr NDCI_low_pass',
          vAxis: {title: 'Cmask and NDCI_gt'+low_pass_th},
          pointSize: 2,
          colors: ['#163bff']});
print(NDCI_low_chart);

var NDCI_mid_chart = ui.Chart.image.seriesByRegion(OUTPUT_images, Sahara_study, 
      ee.Reducer.mean(),'NDCI_mid_pass', 500, 'system:time_start','label')
        .setOptions({
          title: 'Cloudines MODIS_sr NDCI_mid_pass',
          vAxis: {title: 'Cmask and NDCI_gt'+mid_pass_th},
          pointSize: 2,
          colors: ['#fff23b']});
print(NDCI_mid_chart);

var RED_high_chart = ui.Chart.image.seriesByRegion(OUTPUT_images, Sahara_study, 
      ee.Reducer.mean(), 'RED_high_pass', 500, 'system:time_start','label')
        .setOptions({
          title: 'Cloudines MODIS_sr RED high pass',
          vAxis: {title: 'Cmask and RED_sr_gt ('+red_high_pass},
          pointSize: 2,
          colors: ['#ff5371']});
print(RED_high_chart);

var px_count_chart = ui.Chart.image.seriesByRegion(OUTPUT_images, Sahara_study, 
      ee.Reducer.count(), 'RED_high_pass', 500, 'system:time_start','label')
        .setOptions({
          title: 'No. of valid pixels(500m res) in study area',
          vAxis: {title: 'valid pixel px_count'},
          lineWidth: 1,
          pointSize: 2,
          colors: ['#ff5371']
        });
print(px_count_chart);

var MODIS_SR_day = ee.ImageCollection('MODIS/061/MOD09GA')
                                    .filterDate(ee.Date(OUTPUT_Day.get('system:time_start')));
                                    
Map.addLayer(MODIS_SR_day,{'bands':['sur_refl_b01','sur_refl_b04','sur_refl_b03'],'min':100,'max':8000},'rgb');
Map.addLayer(OUTPUT_Day.select('Cflag_state_1km').selfMask(),{palette:'#ffffff'},'C_flag_cloud_state' ,false);
Map.addLayer(OUTPUT_Day.select('NDCI_low_pass').selfMask(),{palette:'#163bff'},'NDCI_low_pass');
Map.addLayer(OUTPUT_Day.select('NDCI_mid_pass').selfMask(),{palette:'#fff23b'},'NDCI_mid_pass',false);
Map.addLayer(OUTPUT_Day.select('RED_high_pass').selfMask(),{palette:'#ff5371'},'RED_high_pass');

Discussion

Year 2001 is standing out with the lowest cloudiness; also the least imagery was used for analytics, surprisingly mayority of missing "dates" is in june when cloudinest is second lowest; meaning that missing dates would lower average cloudiness of year 2001.

In this work cloud state(C-flag) from 1000m state_1km band was not used in analytics. Observations show that when cloudiness of study area is not high;then higher values of cloud state cloudiness compared to Low_pass cloudiness indicate high fogginess. So exporting analytics of this band will be included in further work with latest years added. MODIS is probably the most important instrument for earth observation, due to relatively high spatial resolution and legacy of years.