GEOM184 - Open Source GIS

Practical 3

Google Earth Engine and applications to Harmful Algal Blooms (HABs)

2024/2025

Welcome to the third practical for GEOM184 - Open Source GIS. Today we will start working with Google Earth Engine (GEE) with the applied task to estimate Harmful Algal Blooms (HABs) in Lough Neagh using remote sensing.

Today, we will focus on familiarising with GEE, import remote sensing data, and apply simple func- tions for the estimation of Cyano-HABs indices.

We will use the JavaScript version of the Earth Engine (which is how it was originally developed) . Although this  means  learning  a  new  coding  language,  JavaScript.  is  fairly  intuitive  and  working through this practical will help you familiarise with the basic operations.  Also, remember that if you feel unsure, you can utilise LLMs for support with code editing (for which they normally perform quite well) .

Before starting, you need to make sure that you have access to the  Earth  Engine:  on https: //code.earthengine.google.com/register log in with a Google account (or create it, if you don’t have one) . Then, proceed to Register a Noncommercial or Commercial Cloud Project, select Unpaid usage (as you are using this for education purposes), and select Academia and research from the dropdown menu.  Then select  Create a new Google Cloud Project and continue to summary. Check that everything is working and your account is now ready for use.

Important:   today  we  will  make  use  of  remote  sensing  indices  derived  from  the existing literature.  For some references that you can use for your analysis please refer to the reading list on ELE at https://ele.exeter.ac.uk/course/view.php?id=21417§ion=8.

1    Part A - Getting started and loading remote sensing

In this section, we  prepare our environment and  load  remote  sensing data,  applying filters and extracting several time frames.

1.1    Getting started

First, because Google Earth Engine (from now on GEE for simplicity) is cloud-based, we will work on the online code editor.  Head to https://code.earthengine.google.com/ and log in with your credentials. Whenever possible, save your script by clicking on the Save drop-down menu and choosing Save as. . . .

As you may know from previous courses that you attended, commenting is an essential part when using code-lines.  Try and use proper commenting while working in the code editor  (comments in JavaScript can be made by using //) . At each step we are going to take in this practical, you can run your code to see if your results are displaying correctly and/or error messages appear.

As a first step, we want to include a shapefile representing Lough Neagh in our working environ- ment, which is available from your coursework material.  To add it to GEE, click on the tab Assets on the left-hand side of the code editor and from the dropdown list New select Shape files.  Select your files (you may need to exclude the file ending with  .qmd extension) and then click on Upload. This may take a few seconds, but shortly your new asset should appear.

In your editor window, you can now import your asset and include it within your working environ- ment:

//  Import  Lough  Neagh  polygon

var NI = ee.FeatureCollection( ”projects/ee-YOURNAMEHERE/assets/LoughNeagh ”).geometry();

Be mindful to adjust the pattern above, as it needs to reflect your username and working folder within the code editor.  To visualise  it within our working environment, we can add the following lines:

Map.addLayer(NI, {}, ”LoughuNeagh ”);

Map.centerObject(NI, 10); //  Focus  on  the  region  of  interest

In essence, the command Map .addLayer will add an element to the map, in our case the variable we defined as purpleNI, whilst  ”Lough  Neagh”  it  is  simply the  layer name as it appears on the layer manager of GEE map. The next line is simply a command to centre our map to the area of interest, with a zoom level of 10:  make  it  bigger  (e.g.,  12) and it will zoom in, make it smaller (e.g., 8) and it will zoom out.  If you now run the code, you should be able to see the Lough Neagh polygon overlain to the GEE basemap.

Task 1

Is the Lough Neagh polygon appearing on screen?  If you open the layer manager in the main map, can you change the transparency?

Nevertheless, the particularly complex nature of this shapefile may have detrimental effects on our calculation (so much that it could exceed GEE allocated memory for user usage) .  Therefore, it might be easier to use a bounding box from our shapefile to perform most analysis, and then clip our result to the more complex shapefile. To do so, we add:

//  Compute  bounding  box  (rectangle)  from  Lough  Neagh  polygon

var bbox = NI.bounds();

1.2    Load remote sensing data

We now want to import remote sensing data.  We will work with multi-spectral imagery of Sentinel- 2.  First,  let’s analyse the dataset:  in the search  bar, type Sentinel-2 and when results are loaded select Harmonized Sentinel-2 MSI: MultiSpectral Instrument, Level-2A; this is the Sentinel-2 har- monised datasets for surface reflectance, which is suitable for the type of analysis we will be making. In general, you can look up any type of remote sensing that is available in GEE in this way, so you will know whether you can utilise such catalogue or not.  Next, we can add the full catalogue to our environment:

//  Import  Sentinel-2  level-2A  harmonised  dataset

var S2 = ee.ImageCollection( ”COPERNICUS/S2˙SR˙HARMONIZED ”);

Warning: due to encoding, sometimes copying and pasting directly from this PDF may not show special characters (e.g., the underscore in the above example).  Double check that this appears when you paste your code and edit where necessary.

In this line of code we defined the variable S2 as the collection of Sentinel-2 Level-2A harmonised data.  However, this means that the entire dataset of Sentinel-2 for the whole world is available, and we will need to filter these data in space and time.

The first filtering that we can perform is to restrict our analysis to Lough Neagh.  To do so, we can use several filtering options:

var S2   NI = S2 .filterBounds(bbox)

This function applies a filter that a collection that intersects a geometry, in our case the bounding box of Lough Neagh.  Next, as we are  using multi-spectral data, we want to filter out tiles where cloud cover is excessive.  To do this, we simply attach to the  previous  line of code the following (after starting a new line):

.filterMetadata( ' CLOUDY˙PIXEL˙PERCENTAGE ' ,  ' less˙than ' , 20) //Cloud  cover  filtering   By changing the numerical value, you will be able to change the percentage filter level that we are applying.

2    Defining indices and time filtering

2.1    Define time frames

Here, we will define the time frame in which we want to perform our analysis.  In this practical, we will limit our analysis to the 2020-2024 period and for quarter intervals.

var years = ee.List.sequence(2020, 2024);

var quarters = ee.List.sequence(1, 4);

These two functions only create sequences of integer numbers:  in the first case from 2020 to 2024, and in the second case from 1 to 4 .

Next, for each quarter we need to get a representative value for each pixel. To do so, we will take the median value for each band, for each pixel:

//  Define  a  function  to  calculate median  for  a  given  time  period

var getQuarterlyMedian = function(startDate, endDate) {

var quarterlyImage = ee.Image(S2 NI.filterDate(startDate, endDate).median()).clip(NI); return quarterlyImage;

};

This is a preparatory function that will then be called within a loop.  The function simply estimates the median value for each pixel over each quarter.  Interestingly, notice that we added  .clip(NI): this will make sure that after the analysis, each output will be clipped to the Lough Neagh polygon.

Task 2

Can you think why we would prefer using the median here?  Can you create a similar function for, e.g., average?

2.2    Create a function for environmental indices

We can now start creating our fist function to calculate remote sensing-based indices.  We will start with a simple one:  NDVI. GEE has in-built function to calculate NDVI, simplifying the process.

//  Define  a  function  to  calculate median  NDVI  for  a  given  time  period

var getQuarterlyNDVI = function(image) {

var ndvi = image.normalizedDifference([ ' B8 ' ,  ' B4 ']); return ndvi.clip(NI);

};

What we have done here is creating a new function getQuarterlyNDVI to be executed later; the content within curly brackets { } is the function that will be executed when getQuarterlyNDVI is called; the image .normalizedDifference is used to compute the normalised difference of two variables for the image within GEE to which it is applied. The image .normalizedDifference is an in-built GEE function that effectively computes  x+y/x-y , where (for NDVI) x is B8 near-infrared NIR and y is B4, Red. Once again, please note our clipping at the end of the defined function, as this means that the image is explicitly clipped after all operations have been performed.

We are not really visualising anything here:  rather, we are preparing some key-functions for our next set of analysis.

Task 3

Using the same code structure as above can you generate a function for the  Normalised Difference Cyanobacteria Index (NDCI), e.g., according to Lomeo et al. (2024)?

Of course, we can create custom functions that are not just pre-built into GEE. This could be the case for the Cyanobacteria Index (CI) as defined by e.g., Mishra et al. (2019), which is:

where R , G, and NI R are the red, green and NIR bands, respectively, and λ indicates the middle wavelength for each band. This can be implemented into GEE as follows:

//  Define  constant  values  for wavelengths var l3 = 560; //  Green wavelength  (nm)

var l4 = 665; //  Red wavelength  (nm)

var l8 = 842; //  Near-infrared wavelength  (nm) //  Define  the  function  to  calculate  CI

var calculateCI = function(image) { var CI = image.expression(

' -uB4u+uB3u-u(B3u-uB8)u*u(l4u-ul3)u/u(l8u-ul3) ' , { ' B3 ' : image.select(' B3 '), //  Green  band

' B4 ' : image.select(' B4 '), //  Red  band

' B8 ' : image.select(' B8 '), //  Near-infrared  band ' l3 ' : l3,

' l4 ' : l4,

' l8 ' : l8

});

return CI.rename(' CI ').clip(NI); };

Take a moment to quickly examine the above: apart from defining the wavelength values, we have defined a new function, i.e., CI, using the image .expression that uses a pixel-wise computation.

Task 4

Using a similar approach, can you define the ABDI index by Cao et al. (2021)? This will be necessary to compare different indices in your analysis. Also note the type of bands used by Cao et al. (2021) and their spatial resolution.

3    Part C - Analyse and visualise the map elements

We now need to run our functions that we have created for each quarter and add each iterated value to the map. We will start with NDVI only, then you will add the other indices on your own.

3.1    Create palettes

To visualise colour scales correctly, we need to define palettes. A simple one for NDVI could be:

var paletteNDVI = [ ' red ' ,  ' white ' ,  ' green '];

var visNDVI = {min: -1, max: 1, palette: paletteNDVI};

which is, in essence, a gradient of red to green varying between -1 and 1 (which are the min and max values of NDVI that can be obtained)

Task 5

Can you create similar palettes for the other indices that you have defined?  Bear in mind the potential maximum and minimum values that each index can achieve and design your palettes accordingly.

3.2    Create loop for quarterly calculations

We are now ready to calculate the median pixel-wise NDVI index for each quarter, and to add it to the map for visualisation. We can use the following loop:

years.evaluate(function(yearsList) {

quarters.evaluate(function(quartersList) {

yearsList.map(function(year) {

quartersList.map(function(quarter) {

var startDate = ee.Date.fromYMD(year, (quarter - 1) * 3 + 1, 1);  var endDate = startDate.advance(3,  ' month ').advance(-1,  ' day '); var label =  ' Q '  + quarter +  ' - '  + year;

var image = getQuarterlyMedian(startDate, endDate);

var ndvi = getQuarterlyNDVI(image);

Map.addLayer(ndvi, {min: -1, max: 1, palette: paletteNDVI}, label +  ' -NDVI '); });

});

});

})

Run the code and observe the results appearing on your map.  Under the layer manager, you can select or deselect any quarter of your choice.  You may also get some errors whilst running it.  This may occur if you have set up a too conservative filter (e.g., cloud cover  <5%) and therefore it cannot  be  possible to find any suitable  remote sensing  image for a certain  period, or a  median cannot be established.  In this case to  increase the number of images available you can relax your criteria, if necessary.

Task 6

Why is the majority of the observed values in red (i.e., negative NDVI) for most images? Hint:  search for NDVI ranges and meaning.  Can you find an  instance in which large areas of Lough Neagh had NDVI>0, i.e., there is abundance of green colour being shown? What does this mean in the context of Cyano-HABs?

Now, perhaps you may want to try a different index, to see what information you can obtain.

Task 7

Either by simply changing the name on the loop, or by adding new lines with the indices that you want to test, edit the loop to include, e.g., NDCI, CI, ABDI. Which one seems to give you the right type of information?  Any indication that HABs were occurring at any one time?

4    Conclusions

This is the end of Practical 3 .  By now, you should be confident in using Google  Earth Engine for loading remote sensing data, filtering in space and time, and perform some index analysis.  Next week we will include a few more interactive elements, as well as the ability to produce plots from the images generated in the code editor.

A few points for reflection, considering that this material will be used for your coursework:

      •  Are quarters the right temporal frequency? Can we use something different (e.g., months)?

•  Some indices are suggested, but there are many available in the literature:  can you try and test a few?

•  We limited our analysis between 2020-2024, can we try and increase our time range to, e.g., 2017-2024?

• Although we will work on this next time, can you start thinking of some of the quantitative analysis that you can show through these initial results?