代写GEOM184 - Open Source GIS Practical 4 2024/2025调试数据库编程

日期：2025-03-15 05:36

GEOM184 - Open Source GIS

Practical 4

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

Welcome to the fourth practical for GEOM184 - Open Source GIS. Today we will continue working with Google Earth Engine (GEE) and the Harmful Algal Blooms (HABs) problem.

Today, we will mostly focus on including additional elements to our GEE visualisation and how to produce quantitative plots of the indices that we have defined last week.  We will also  learn how to export the data produced, and to publish our analysis as an app for public use (this is not a requirement for your assessment, but it is a really helpful tool) .  Remember that if you feel unsure, you can utilise LLMs for support with code editing (for which they normally perform. reasonably well) .

Important:  we will use the code that we have generated last week, so if you are not up to date with it, this is the time to complete it!  Also, please note that due to formatting issues, some of the code lines may need little tweaks, e.g., missing underscores.

1 Part A - Creating a split panel

In this section, we will learn on how to create a split panel:  this will give us a chance to visualise side-by-side two indices of our choice (or more, by switching between layers) .  In this way, we will get a true understanding of the occurrence of Cyano-HABs in Lough Neagh.  We will keep most of the code that we produced last week, but we will need to make some changes within the loop and add a few lines just before the loop itself.

1.1 Remove some of the old lines

We need to remove a few lines to get the split panel to work. You can remove (or, if you prefer, comment out by using //) the following:

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

Map.centerObject(NI, 10);

1.2 Create two panels

Now, we want to create two map instances, i.e., we are creating two separate map environments. To do this, include the following lines of code within this subsection before the loop, (i.e., the part of the code beginning with years.evaluate):

//  Create  two map  instances

var map1 = ui.Map();

var map2 = ui.Map();

This code creates two map instances within the GEE environment, and each map environment (i.e., map1 and map2) can have elements added to it independently.  The ui.Map() is a function of GEE that serves the creation of an interactive map.  Each of the newly defined maps will work independently, unless explicitly linked (which we will do soon) .

Next step, we centre both our new maps to Lough Neagh, so we add:

//  Set  center  and  zoom  level

map1 .centerObject(NI, 10);

map2 .centerObject(NI, 10);

You can notice that there is very little difference between these lines and the one we deleted earlier, except for explicitly mentioning map1 and map2.  Now, the next step is to add the maps to the user interface, where we assign equal width to both maps (and 100% height) .  To do so  include the following lines:

//  Add  the maps  to  the  user  interface

var panel1 = ui.Panel({

widgets: [map1],

style.: {width:  ' 50% ' ,  height:  ' 100% '}

});

var panel2 = ui.Panel({

widgets: [map2],

style.: {width:  ' 50% ' ,  height:  ' 100% '}

});

The ui .Panel function  requires two arguments:  widgets  (where we select the  map  for each panel), and style where we can change the appearance (and in our case decide width and height of the panels) .

Next, we add the panels to the split panel:

//  Add  the  panels  to  a  split  panel

var splitPanel = ui.SplitPanel({

firstPanel: panel1,

secondPanel: panel2,

wipe: true

});

And we can now add the split panel to the user interface:

//  Add  the  split  panel  to  the  user  interface

ui.root.clear();

ui.root.add(splitPanel);

All of the above is great, but now we need to make sure that the two panels are linked to each other, so that if any operation is performed (zooming in and out, panning, etc.)  will be observed in both panels. So we add a map linker instance:

//  Create  a  Map  Linker  instance

var mapLinker = ui.Map.Linker([map1, map2]);

1.3    Adapt the loop

All is ready now to visualise a split panel with the indices of your choice.  However, to do so we need to slightly adapt our loop.  For example, last week we used the following for NDVI:

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

Assuming we want on one side NDVI and on the other side another index (e.g., ABDI), we can remove (or comment out with //) the above code and add instead:

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

map2 .addLayer(ABDI, {min: -1, max: 1, palette: paletteABDI}, label +  ' ˙ABDI ');

Remember that any index you are adding within the loop needs to be defined in the same way as var  NDVI  =  getQuarterlyNDVI(image); and this needs to occur before you add the new layers to either map1 or map2.  Now, try to run your code and see how it looks with the results.

Task 1

Can you move the slider left and right?  Are the two maps joined together in the same place and remained joined even when zooming in and out?  Can you tick and untick layers for quarters of your choice for each map?  (hint:  you  may have to move the slider all the way to the right to open the layer manager for the left map)

1.4    Explore visualisation of different indices

Now that you have two maps to visualise your results, we can use them to make a comparative analysis between quarters. You may want to focus on Summer 2023 as a prime example oh Cyano- HABs observations  (especially due to the exceptionality of the event) .   One way to do it is to analyse indices pair-wise.  It  is of course possible to add more indices to your maps.  This is up to you in what way you want to perform your analysis.  For example, you may want to add NDVI and NDCI for map1 and CI and ABDI to map2 - to do so, just follow the same steps as above and make sure to add layers to the map of your choice.

Task 2

Can you see any sign of Cyano-HABs in the map? Can you refer back to the interpretation of these indices in the literature (e.g., NDVI>0 indicates vegetation, and the nearer it gets to 1, the healthier the vegetation, and so if you would expect vegetation in the middle of a lake)?

2 Part B - Statistical analysis of indices

To have a full understanding and a quantitative assessment of Cyano-HABs in Lough Neagh, we may need to consider a (relatively basic) statistical analysis of data that can quantify Cyano-HABs occurrence, for which we  may want to create  plots directly within  GEE. Of course,  you could download individually the maps that you have created and analyse each image within QGIS, but this would not be a good use of your time and of computational resources.  To this end, we need to edit our code in two steps: just before the loop, and within the loop.

2.1 Preparation for the loop

This step considers that before the loop we need to introduce some empty lists that will be populated by the numerical values generated at each iteration. Assuming we are using NDVI and ABDI, the first part of the code to include is:

//  Lists  to  hold  NDVI  and  ABDI  stats

var ndviValues = [];

var abdiValues = [];

This simply adds two empty lists where we will store the information for both NDVI and ABDI.

2.2    Amending the loop

We can now make a few changes within the loop to produce an output of our statistics in the form of a plot that includes mean and set percentile values for NDVI (and ABDI), for each quarter.

In the most nested part of the loop (just after defining variables NDVI and ADBI) we add:

//  Compute  statistics  for  NDVI

var ndviStats = ndvi.reduceRegion({

reducer: ee.Reducer.mean()

.combine(ee.Reducer.percentile([5, 95]),  ' ' , true),

geometry: NI,

scale: 10

});

The newly defined function ndviStats is the one that contains all statistical values that we are using for our analysis.  The  method utilised for this function is the so-called reduceRegion, and we then define a reducer: this is an identification for the operation that we want to perform.  For example, the first argument is ee.Reducer.mean(), which simply means that we are calculating the mean value of NDVI for the whole of Lough Neagh for each quarter (don’t forget that this operation lies within the most nested part of the loop, meaning that each operation is repeated for each loop iteration, i.e., for every quarter) .  Then, ee.Reducer.percentile is for us to define the NDVI value for that quarter across the whole of Lough Neagh that is, respectively, 5% and 95% of the overall observations. Other helpful percentiles that we could potentially use are 25% and 75% (you may known as the interquartiles) - you can try one and then the other to see what makes the most statistical sense.

Task 4

Can you add the function computing statistical indicators for ABDI (and/or any other index) to this?

Next, is to add the statistical values we have calculated to our empty lists (remember we defined them earlier?):

//  Add  statistics  to  the  lists

ndviValues.push(ee.Feature(null , {

label: label,

Mean: ndviStats.get(' nd˙mean '),

p5: ndviStats.get(' nd˙p5 '),

p95: ndviStats.get(' nd˙p95 ')

}));

Here, be mindful of the names within the function arguments: values on the left hand-side are the newly defined values, whist on the righ-hand side are being pulled from the ndviStats variable that we defined earlier.  If you use different percentiles, you will need to edit this. Important:  please be mindful of the nomenclature used to call the variables stored by ndviStats:  we use the name nd because NDVI was generated using the in-built function of normalised difference within GEE. For any other index that does not use a normalised difference (e.g., ABDI), you will need to use the name that you have renamed when defining calculateABDI (or similar), for which we may have used return  ABDI.rename( ' ABDI ' ) .clip(NI); - therefore we need to replace nd with whatever name we included within ” for calculateABDI.

Task 5

Can you repeat the same operation for ABDI (and any other index you have included) to this?  Please bear in mind the type of function you used (i.e., a normalised difference or a custom function) and the name you may have assigned (see above) .

One final step and we are ready to add mean NDVI and ABDI values as a plot box.  To do this, we need to jump two levels of the loop, meaning that the following lines of code need to be placed within the quarters .evaluate level of the loop, i.e., after the whole set of quartersList .map and yearsList .map iterations have been completed within the loop:

//  Create  NDVI  Chart

var ndviChart = ui.Chart.feature.byFeature(ee.FeatureCollection(ndviValues),  ' label ' ,

[ ' Mean ' ,  ' p5 ' ,  ' p95 '])

.setChartType(' LineChart ')

.setOptions({

title:  ' NDVIuoverutime ' ,

hAxis: {title:  ' Quarter '},

vAxis: {title:  ' NDVIuValue '},

lineWidth: 1,

pointSize: 4

});

Note that in the first line of this code within the square brackets we included the statistical quantities we defined earlier.  You will  need to edit this  if you change the percentiles.  The meaning of the elements inside the other lines of code is relatively simple:  we are plotting a chart using features of a collection, and we are setting options within the chart - feel free to edit any of these to suit your design plan.

Task 6

Can you add ABDI (and any other index) to this?  Would you try to  change some of the visualisation options?

Finally, just after these new lines, we can finalise the creation of a panel within GEE environment:

//  Create UI  panels  for  the  charts

var ndviChartPanel = ui.Panel({

widgets: [ndviChart],

style.: {width:  ' 500px ' , height:  ' 200px ' , position:  ' bottom-left '}

});

//  Add  the  chart  panels  to  the map  panels

map1 .add(ndviChartPanel);

Task 7

Can you add ABDI (and any other index) to this?  Bear in mind that if you followed this approach so far of adding  NDVI to  panel  1  and  ABDI to  panel  2,  ABDI will  be  on the right-hand side panel, so you may want to change the style options!

If you run your code now, you should be able to visualise the plots appearing on screen, indicating the mean and any other percentile that you may have used for your analysis.  If you want more details on these plots, you can click on the icon in the top-right of each plot box, which will open a new window with your plot.  In this  new window you can save the plot as it is, or download the data.

2.3    Analysis and anomaly detection

As we want to understand how severe Cyano-HABs have been for Lough Neagh, the graphs can already tell us some interesting stories.

Task 8

Can you visually see any pattern with the indices?  Any  indication that something unusual occurred in 2023?  Make sure you comment or add information about this.

An even better approach than just visually observing trends is to use a simple Z-score approach. This is given by the ratio:

Where Z is an individual value (e.g., the mean NDVI for a quarter of choice), µ is the average of all these values (e.g., the average of all mean NDVI across all quarters) and σ is its standard deviation. Typically, any value >1 or <-1 indicate that the data point is considered unusual or far from the average value (e.g., greater than at least one standard deviation); the grater the absolute value of the z-score, the more unusual the data point is.

Task 9

By using the Z-score, can you see immediately any event that seemed anomalous?  Make sure you are performing your analysis on the mean rather than the percentiles. You can use a spreadsheet to do this calculation and you don’t need Google Earth Engine.

3    Part C -  Export your data and publish your app

3.1 Exporting results as raster

Finally, we have identified a period of interest, indices that well represent our analysis and we want to produce a map for our report or further analysis.  To do this, we need to add a few more lines at the end of our code (i.e., after the loop):

//  Define  the  date  range  for  Q3-2023  (1st  July  to  30th  September  2023)

var startDate = ee.Date.fromYMD(2023, 7, 1);

var endDate = ee.Date.fromYMD(2023, 9, 30);

//  Get  the median  image  for  Q3-2023

var image = getQuarterlyMedian(startDate, endDate);

//  Calculate  NDVI  for  the  image

var ndvi = getQuarterlyNDVI(image);

//  Export  the  NDVI  image  to  Google  Drive

Export.image.toDrive({

image: ndvi, //  The  NDVI  image  to  export

description:  ' NDVI˙Q3˙2023 ' , //  The  name  for  the  exported  file

scale: 10, //  The  resolution  (in  m)

region: NI, //  The  region  of  interest  (Lough  Neagh)

fileFormat:  ' GeoTIFF '  //  Export  as  GeoTIFF });

We are effectively repeating some of the operations we have done before, because this time we are targeting the period of time we want to perform our export.  Do this for any index you want to use/display and for any relevant period where you believe Cyano-HABs occurred.  You may need to adjust the resolution if the resulting file is too large.

Once the code is run, you will need to action this within your  Tasks environment:  click on  the ribbon and then click Run for the TIF file that you have requested.  It may take some time, but once completed you can download your file in your Google drive and use it within QGIS for inclusion within a layout.

3.2 Publish your app

This is not necessary for this module, but in the future you may want to unlock the potential of GEE and publish your work and share it widely (and, also, if you wish you can include a link within your coursework, although this will not be assessed) .  To create a new app, click on Apps, and when the new window opens, click on New app.  From there, follow the prompts (e.g., you may want for your app to be editable by yourself only) until the app is published and you have a web link to it. And that’s it!

4 Conclusions

This is the end of Practical 4 .  By now, you should be confident in using Google  Earth Engine for plotting data and analysing trends in time and space, as well as identifying anomalies, such as for the Cyano-HABs case we have been working on. See you next time!