Summary and Best Practice

To wrap the first chapter up, we will have a look at an adjusted, summarized script from the snippets of code as well as a new example on how to download Landsat 8 data of 8 years for a single administrative unit of your choice.

For the script we already know, we will look at two different versions: one without any styling techniques and comments, and one that follows basic ‘best-practice’-rules to see the huge difference in quality for your code.

We will then save the scripts to the “Scripts” tab to always have them available for future reference.

Messy Code

This is an example of how your code should not look like.

Shortcomings:

No comments to reconstruct what is being done
No appropriate spacings for better orientation
No names assigned to prints and maplayers, making it hard to tell them apart
No semicolons to mark finished coherent code

var extent_lebanon = ee.FeatureCollection("FAO/GAUL/2015/level0")
                  .filter(ee.Filter.eq('ADM0_NAME', 'Lebanon'))

var s2a = ee.ImageCollection('COPERNICUS/S2_SR')
                  .filterBounds(extent_lebanon)
                  .filterDate('2020-01-01', '2020-12-31')
                  .select('B1','B2','B3','B4','B5','B6','B7','B8','B8A','B9','B11','B12')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10))
print(s2a)

var visParams = {'min': 400,'max': [4000,3000,3000],   'bands':'B8,B4,B3'}
Map.addLayer(s2a, visParams ,'Sentinel2', false)
Map.centerObject(s2a,9)

var s2a_image = s2a.mosaic()
print(s2a_image)

var s2a_clip = s2a_image.clip(extent_lebanon)
print(s2a_clip)
Map.addLayer(s2a_clip, visParams)

var s2_mean = s2a.reduce(ee.Reducer.mean())
              .clip(extent_lebanon)
print(s2_mean)

var visParams_mean = {'min': 400,'max': [4000,3000,3000],   'bands':'B8_mean,B4_mean,B3_mean'}
Map.addLayer(s2_mean, visParams_mean, 'IC 2020 reduce mean', false)

var s2_minMax = s2a.reduce(ee.Reducer.minMax())
              .clip(extent_lebanon)
print(s2_minMax)

var visParams_min = {'min': 0,'max': [2000,2000,2000],   'bands':'B8_min,B4_min,B3_min'}
var visParams_max = {'min': 1500,'max': [8000,9000,9000], 'bands':'B8_max,B4_max,B3_max'}
Map.addLayer(s2_minMax, visParams_min, 'IC 2020 reduce min', false)
Map.addLayer(s2_minMax, visParams_max, 'IC 2020 reduce max', false)

var chart = ui.Chart.image
  .series({
    imageCollection: s2a.select('B3','B4','B8'),
    region: extent_lebanon,
    reducer: ee.Reducer.mean(),
    scale: 200
     })
  .setOptions({
          title: 'Mean Surface Reflectance Value by Date for Lebanon',
          hAxis: {title: 'Date', titleTextStyle: {italic: false, bold: true}},
          vAxis: {title: 'Surface Reflectance',titleTextStyle: {italic: false, bold: true}},
  })
print(chart)

var doychart = ui.Chart.image
  .doySeries({
    imageCollection: s2a.select('B3','B4','B8'),
    region: extent_lebanon,
    regionReducer: ee.Reducer.mean(),
    scale: 200,
    yearReducer: ee.Reducer.mean(),
    startDay: 1,
    endDay: 365
  })
  .setOptions({
          title: 'Mean Surface Reflectance Value by Date for Lebanon',
          hAxis: {title: 'Day of year', titleTextStyle: {italic: false, bold: true}},
          vAxis: {title: 'Surface Reflectance',titleTextStyle: {italic: false, bold: true}},
  })
print(doychart)

Clean Code

This is an example of how your code should look like.

As you can see, the final results of both scripts are the same. However, for the commented script, it is much easier to comprehend the script and retrace single arguments and outputs. This can be essentially useful not only for others reading your script, but also for yourself to easily get back into it when you have not used it for a while.

//Step 1: Access your boundary-defining geometry
var extent_lebanon = ee.FeatureCollection("FAO/GAUL/2015/level0")
                  .filter(ee.Filter.eq('ADM0_NAME', 'Lebanon')); //filter for entry that equals the UN country name 'Lebanon'



//Step 2: Access the Sentinel-2 Level-2A data and filter it for all the the images of the year 2020 that lie within the geometries boundaries. Keep only the relevant bands and filter for cloud coverage.
var s2a = ee.ImageCollection('COPERNICUS/S2_SR')
                  .filterBounds(extent_lebanon)
                  .filterDate('2020-01-01', '2020-12-31')
                  .select('B1','B2','B3','B4','B5','B6','B7','B8','B8A','B9','B11','B12')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10));

//Print your ImageCollection to your console tab to inspect it
print(s2a, 'Image Collection 2020');



//Step 3: Visualizing Data
//Define the visualizing parameters; for this example, we will only define the visualized bands and the min/max values
var visParams = {'min': 400,'max': [4000,3000,3000],   'bands':'B8,B4,B3'}; //bands 8, 4 and 3 are to be visualized; minimum value is 400 for all three bands, maximum value is 4000 for band 8, 3000 for band 4 and 3

Map.addLayer(s2a, visParams ,'Sentinel2', false);   // Visualize the ImageCollection in the Map panel and assign a name to it; Hide the layer by default, as it takes long to render
Map.centerObject(s2a,9);              // Center the map panel on the ImageCollection; zoom level is 9



//Step 4: Further process your ImageCollection
//Turn the ImageCollection to a single layerstack; this is essential for being able to clip it
var s2a_image = s2a.mosaic();
print(s2a_image, 'Layerstack Image 2020 Lebanon, < 10% Cloud Cover');

//Clip the layerstack to the extent of the geometry
var s2a_clip = s2a_image.clip(extent_lebanon);
print(s2a_clip, 'Layerstack Image 2020 Lebanon clipped, < 10% Cloud Cover');

//Add the result to the map
Map.addLayer(s2a_clip, visParams ,'Sentinel-2 Level-2A 2020 clipped');



//Step 5: Use Reducers to have a look at image statistics via the print tab or the Inspector
var s2_mean = s2a.reduce(ee.Reducer.mean())
              .clip(extent_lebanon);        //The result can be clipped, as it is an Image!
print(s2_mean, 'IC 2020 reduce mean');

//the band names have changed, so we will have to define a new visParams
var visParams_mean = {'min': 400,'max': [4000,3000,3000],   'bands':'B8_mean,B4_mean,B3_mean'};
Map.addLayer(s2_mean, visParams_mean, 'IC 2020 reduce mean', false);

//Same for min max values. Take note that this will produce double the bands, one for min and one for max values
var s2_minMax = s2a.reduce(ee.Reducer.minMax())
              .clip(extent_lebanon);
print(s2_minMax, 'IC 2020 reduce minMax');

var visParams_min = {'min': 0,'max': [2000,2000,2000],   'bands':'B8_min,B4_min,B3_min'};
var visParams_max = {'min': 1500,'max': [8000,9000,9000], 'bands':'B8_max,B4_max,B3_max'};
Map.addLayer(s2_minMax, visParams_min, 'IC 2020 reduce min', false);
Map.addLayer(s2_minMax, visParams_max, 'IC 2020 reduce max', false);



//Step 6: Have a look at various charts for your data
//Add a time series chart of the mean surface reflectance values of bands 3, 4 and 8
var chart = ui.Chart.image
  .series({
    imageCollection: s2a.select('B3','B4','B8'),
    region: extent_lebanon,
    reducer: ee.Reducer.mean(),
    scale: 200
     })
  .setOptions({
          title: 'Mean Surface Reflectance Value by Date for Lebanon',
          hAxis: {title: 'Date', titleTextStyle: {italic: false, bold: true}},
          vAxis: {title: 'Surface Reflectance',titleTextStyle: {italic: false, bold: true}},
  });

print(chart);

//Add a day-of-year-chart of the mean surface reflectance values of bands 3, 4 and 8
var doychart = ui.Chart.image
  .doySeries({
    imageCollection: s2a.select('B3','B4','B8'),
    region: extent_lebanon,
    regionReducer: ee.Reducer.mean(),
    scale: 200,
    yearReducer: ee.Reducer.mean(),
    startDay: 1,
    endDay: 365
  })
  .setOptions({
          title: 'Mean Surface Reflectance Value by Date for Lebanon',
          hAxis: {title: 'Day of year', titleTextStyle: {italic: false, bold: true}},
          vAxis: {title: 'Surface Reflectance',titleTextStyle: {italic: false, bold: true}},
  });

print(doychart);

Example for Landsat 8

Lets have a look at the same concept of script for Landsat 8 data. Please run the code in a separate tab.

To spice things up, we will look at an alternative way to explore research areas and observe a longer time series.

//Step 1: Access your boundary-defining geometry
//load the GAUL-Level2 layer to see administrative units. choose one and place a point geometry in it. This code should be commented as soon as you chose one!
//var gaul_l2 = ee.FeatureCollection("FAO/GAUL/2015/level2");
//Map.addLayer(gaul_l2);

//The point geometry will automatically be added to your Imports at the top of your code
//To make this reconstructable, this is the point that was used for this tutorial. If you chose your own, please skip the following line.
var geometry = /* color: #d63000 */ee.Geometry.Point([38.931862676219545, 35.80741370345301]);



//Filter the GAUL-Level2 layer to only contain the desired feature you placed your point geometry in.
var au_layer = ee.FeatureCollection('FAO/GAUL/2015/level2')
                    .filterBounds(geometry);
print(au_layer, 'Administrative Unit');
Map.addLayer(au_layer, {},'Administrative Unit');
/*
We can now see detailed information on the Administrative Unit we selected by browsing 'features' -> 0 -> properties, for example:
ADM0_CODE: 238
ADM0_NAME: Syrian Arab Republic
ADM2_NAME: Raqqa
*/



//Step 2: Choose a Landsat 8 Catalog for your needs and filter it. In this example, we use Tier 1 Surface Reflectance data. Make sure to check the documentation to learn more about its bands and image properties!
var ls8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
        .filterBounds(au_layer)
        .filterDate('2013-04-1', '2021-04-1')
        .filter(ee.Filter.lt('CLOUD_COVER', 10));
print(ls8, 'Image Collection 2013-2021 <10% clouds'); //we receive 439 elements



//Step 3: Define the visualizing parameters
var visParams = {'min': 400,'max': [4000,3000,3000],   'bands':'B5,B4,B3'}; //bands 5, 4 and 3 are to be visualized; minimum value is 400 for all three bands, maximum value is 4000 for band 5, 3000 for band 4 and 3

Map.addLayer(ls8, visParams ,'Landsat 8', false);   // Visualize the ImageCollection in the Map panel and assign a name to it; Hide the layer by default, as it takes long to render
Map.centerObject(ls8,10);              // Center the map panel on the ImageCollection; zoom level is 10



//Step 4: Further process your ImageCollection
//Turn the ImageCollection to a single layerstack; this is essential for being able to clip it
var ls8_image = ls8.mosaic();
print(ls8_image, 'Layerstack Image 2013-2021, < 10% Cloud Cover');

//Clip the layerstack to the extent of the geometry
var ls8_clip = ls8_image.clip(au_layer);
print(ls8_clip, 'Layerstack Image 2013-2021 clipped, < 10% Cloud Cover');

//Add the result to the map
Map.addLayer(ls8_clip, visParams ,'Landsat8 T1 SR 2013-2021 clipped');



//Step 5: Use Reducers to have a look at image statistics via the print tab or the Inspector
var ls8_mean = ls8.reduce(ee.Reducer.mean())
              .clip(au_layer);        //The result can be clipped, as it is an Image!
print(ls8_mean, '2013-2021 reduce mean');

//the band names have changed, so we will have to define a new visParams
var visParams_mean = {'min': 400,'max': [5000,5000,5000],   'bands':'B5_mean,B4_mean,B3_mean'};
Map.addLayer(ls8_mean, visParams_mean, '2013-2021 reduce mean', false);

//Same for min max values. Take note that this will produce double the bands, one for min and one for max values
var ls8_minMax = ls8.reduce(ee.Reducer.minMax())
              .clip(au_layer);
print(ls8_minMax, '2013-2021 reduce minMax');

var visParams_min = {'min': 0,'max': [2000,2000,2000],   'bands':'B5_min,B4_min,B3_min'};
var visParams_max = {'min': 1500,'max': [8000,9000,9000], 'bands':'B5_max,B4_max,B3_max'};
Map.addLayer(ls8_minMax, visParams_min, '2013-2021 reduce min', false);
Map.addLayer(ls8_minMax, visParams_max, '2013-2021 reduce max', false);



//Step 6: Have a look at various charts for your data
//Add a time series chart of the mean surface reflectance values of bands 3, 4 and 5
var chart = ui.Chart.image
  .series({
    imageCollection: ls8.select('B3','B4','B5'),
    region: au_layer,
    reducer: ee.Reducer.mean(),
    scale: 200
     })
  .setOptions({
          title: 'Mean Surface Reflectance Value by Date',
          hAxis: {title: 'Date', titleTextStyle: {italic: false, bold: true}},
          vAxis: {title: 'Surface Reflectance',titleTextStyle: {italic: false, bold: true}},
  });

print(chart);

//Add a day-of-year-chart of the mean surface reflectance values of bands 3, 4 and 5
var doychart = ui.Chart.image
  .doySeries({
    imageCollection: ls8.select('B3','B4','B5'),
    region: au_layer,
    regionReducer: ee.Reducer.mean(),
    scale: 200,
    yearReducer: ee.Reducer.mean(),
    startDay: 1,
    endDay: 365
  })
  .setOptions({
          title: 'Mean Surface Reflectance Value by Date',
          hAxis: {title: 'Day of year', titleTextStyle: {italic: false, bold: true}},
          vAxis: {title: 'Surface Reflectance',titleTextStyle: {italic: false, bold: true}},
  });

print(doychart);

We can now save our finished scripts to our script repository.

To do so, simply click on Save in the Code Editor and choose a path and a file name. Like this, we always have them available for future reference on how to access, filter, clip and visualize S-2 Level-2A or Landsat 8 Tier 1 SR data for our research area. All we have to change is adjusting the research area as well as choosing the respective timespan.

It might be useful to build up a folder structure similar to organizing data on your computer, as scripts can easily add up. When you have separate folders for specific sensor- or applications-code, you will have an easier time finding, managing and organizing your scripts.

Department of Earth Sciences

GEO-IT - The Technology of Data Acquisition for Sustainable Development and Crisis Management

Summary and Best Practice

Messy Code

Clean Code

Example for Landsat 8