11 Cleaning up the Imagery
One of the primary problems with using satellite imagery is the location of our sensor, that’s recording the data, is in space! That, in itself, isn’t much of a problem. However, the fact that it’s above the clouds is a problem. Passive sensors, that record the sun’s reflectance, generally have a very difficult time seeing through the clouds. Yes, there are some active sensors that can accomplish this feat, but that’s out of the scope of this tutorial. Now back to the issue at hand, we need to get rid of the clouds.
Historically, it was pretty much sheer luck if you could find an image without cloud cover. Clouds are a major part of our water cycle and particularly important in moving water from one place to another. Modern data providers, like the European Space Agency (ESA), have filters in their search engines to help identify images with fewer clouds. While this is effective, GEE utilizes algorithms to find reflectance data without clouds…when possible. Try finding a cloud free month over the United Kingdom. It may not happen. Sometimes, clouds happen…relax.
In this section, we’ll take a look at our two major sensors – Sentinel and Landsat. Although these represent a series of sensors, we’ll use some of the more popular vehicles.
Sentinel 2
You’ve already been using the ESA’s Sentinel 2 data in the previous section, so you already know a couple of things:
- The location of Sentinel 2 data in the EE Data Catalog is accessed through the following function ee.ImageryCollection(‘COPERNICUS/S2_SR_HARMONIZED).
- To display the data, you have two options:
- Pointing to an exact image in the catalog (you haven’t done this)
- Filtering by a series of dates and showing the mean reflectance using the .filterDate() and .mean() functions.
- You must define the bands to display through a visual parameters object using the following format:
- var unique_visparams = {bands:[‘B4’, ’B3’ ,B2’], min: 0, max: 4000};
- The Map.addLayer() expects the following format:
- Map.addLayer(data_variable_name, visparams, ‘name in display’);
Hopefully that was a useful refresher and a guide you can use in the future. Now, the smart people over at GEE helped develop a function to remove as many clouds as possible. We haven’t seen a function yet, but it’s an incredibly useful way of defining a tool that can be applied at any point in the script. Let’s take a look at a simple example of defining a function in JavaScript:
Example 3.1
- function addstuff(thing1,thing2) {
- print(thing1+thing2);
- }
- addstuff(10,12);
You can see in the example above, how we create a function that adds stuff together. We require two parameters: thing1 and thing2, which the user defines when they call the function. Although we could have just said print(10+15), this function could be a complex series of instructions that could be applied over a dataset as many times as required, saving us a boatload of time. Here’s a slightly more complicated use of the same function:
Examples 3.2
- function addstuff(thing1,thing2) {
- print(thing1+thing2);
- }
- var lots_of_numbers = [2, 4, 8, 10, 14, 16, 20, 22, 26, 28];
- var total_numbers = lots_of_numbers.length;
- for(var i=0; i < total_numbers; i++) {
- addstuff(lots_of_numbers[i],lots_of_numbers[i]);
- }
By using the function you save having to write print(lots_of_numbers[n], lots_of_numbers[n]) ten times! Already a huge time savings. Yeah, I threw a “for loop” at you. Don’t worry…relax. You can handle this. We’ll save that explanation for another time, but suffice it to say we’re just using that loop to iterate our way through that list of numbers.
Now how does this all relate to removal of clouds in Sentinel 2 imagery? Because we remove the clouds through a function. See, Sentinel 2 data comes with a ‘QA’ Band that identifies clouds and cirrus with an identifier. In this case, those pixels get a number 10 (bit) and 11 (bit) respectively. Pixels that don’t have clouds get a (bit) value of 0. The function plays a little game to move the cloudy pixels to a (bit) value of 1 and clear conditions get a value of 0. Here’s what that function looks like (directly from the EE Data Catalog Example):
Example 3.3
- /**
- * Function to mask clouds using the Sentinel-2 QA band
- * @param {ee.Image} image Sentinel-2 image
- * @return {ee.Image} cloud masked Sentinel-2 image
- */
- function maskS2clouds(image) {
- var qa = image.select(‘QA60’);
- // Bits 10 and 11 are clouds and cirrus, respectively.
- var cloudBitMask = 1 << 10;
- var cirrusBitMask = 1 << 11;
- // Both flags should be set to zero, indicating clear conditions.
- var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
- .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
- return image.updateMask(mask).divide(10000);
- }
Let’s apply this to the example to one of our Part II Exercises, where we look at the Swiss Alps. This is, so far, the largest code you’ve seen so far. Just take it line by line:
Example: Sentinel 2 Cloud Removal Function
- //Create the maskS2clouds function
- function maskS2clouds(image) {
- var qa = image.select(‘QA60’);
- // Bits 10 and 11 are clouds and cirrus, respectively.
- var cloudBitMask = 1 << 10;
- var cirrusBitMask = 1 << 11;
- // Both flags should be set to zero, indicating clear conditions.
- var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
- .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
- return image.updateMask(mask).divide(10000);
- }
- //Bring in the correct collection and filter by the desired date.
- var imagery_data = ee.ImageCollection(‘COPERNICUS/S2_SR_HARMONIZED’)
- .filterDate(‘2023-05-01′,’2023-05-30’)
- .filter(ee.Filter.lt(‘CLOUDY_PIXEL_PERCENTAGE’,70))
- .map(maskS2clouds);
- //Take the mean reflectance values for the imagery taken during the month
- var mean_image = imagery_data.mean();
- //Establish visual parameters for the imagery to display standard, visible imagery
- var RGB_params = {
- bands:[‘B4′,’B3′,’B2’],
- min: 0.0,
- max: 0.4
- };
- //Center the map somewhere over the middle of the alps. Zoom 7.5x.
- Map.setCenter(9.664, 46.4458, 7.5);
- //Add the layers to map
- Map.addLayer(mean_image, RGB_params, ‘RGB Imagery’);
Take a breath. That’s a lot to chew. The first 15 lines just define the maskS2clouds function. You don’t need to memorize this. Just copy and paste it into your script each time you need to remove clouds from Sentinel2 imagery. A couple of things you need to remember:
- You’ll need to define the maximum amount of pixels that can be clouds with the filter less than (lt) function in your imagery pull. Normally I start at 20% cloud coverage but the alps had a lot of clouds in May of 2023. I had to keep increasing my percentage to 70% to get a full mosaic:
- .filter(ee.Filter.lt(‘CLOUDY_PIXEL_PERCENTAGE’, 70))
- You need to pass the mosaic through the maskS2clouds() function, when you pull the data:
- .map(maskS2clouds)
- Notice that the min, max has been reduced to 0.0 and 0.4 respectively? It’s really a pain dealing with the big numbers. The maskS2clouds() function divides all the values by 10,000. We’ll do that from here out.
That’s it. Go ahead and add the glacier dataset back in. I’m going to change the color of this dataset to blue to make it easier to see:
Example 3.4
- //Load the GLIMS dataset for glaciers and filter for the appropriate geographic region
- var dataset = ee.FeatureCollection(‘GLIMS/current’).filter(‘geog_area == “Swiss Alps”‘);
- //Establish visual parameters for the glacier data.
- var alps_visparams = {
- color:’blue’
- };
- //Center the map somewhere over the middle of the alps. Zoom 7.5x.
- Map.setCenter(9.664, 46.4458, 7.5);
- function maskS2clouds(image) {
- var qa = image.select(‘QA60’);
- // Bits 10 and 11 are clouds and cirrus, respectively.
- var cloudBitMask = 1 << 10;
- var cirrusBitMask = 1 << 11;
- // Both flags should be set to zero, indicating clear conditions.
- var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
- .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
- return image.updateMask(mask).divide(10000);
- }
- //Bring in the correct collection and filter by the desired date.
- var imagery_data = ee.ImageCollection(‘COPERNICUS/S2_SR_HARMONIZED’)
- .filterDate(‘2023-05-01′,’2023-05-30’)
- .filter(ee.Filter.lt(‘CLOUDY_PIXEL_PERCENTAGE’,70))
- .map(maskS2clouds);
- //Take the mean reflectance values for the imagery taken during the month
- var mean_image = imagery_data.mean();
- //Establish visual parameters for the imagery to display standard, visible imagery
- var RGB_params = {
- bands:[‘B4′,’B3′,’B2’],
- min: 0.0,
- max: 0.4
- };
- //Add the layers to map
- Map.addLayer(mean_image, RGB_params, ‘RGB Imagery’);
- Map.addLayer(dataset, alps_visparams, ‘Swiss Alps’);
Now apply what you’ve learned here to question 1 of Exercises 2.1 on the Badger Fire for July of 2023. See how low of a percentage you can go before your mosaic is missing pieces. After you’re done with this, we’ll take a similar peak at Landsat 8 imagery.
Landsat 8
Landsat is to NASA as Sentinel is to the ESA. Landsat is a joint operation between the United States Geological Survey (USGS) and NASA and has been observing the Earth since 1972 (https://developers.google.com/earth-engine/datasets/catalog/landsat). Currently, GEE has the 9th generation of satellite data available. As this data gets fully ingested, we’ll focus on Landsat 8.
To bring in the Landsat 8 Surface Reflectance Imagery, it’s pretty much the same routine as Sentinel. A couple of weird things:
- The range of values is strange: 0 to ~15,000. I actually don’t know what the upper limit is, but it’s in that ballpark. Because of this we apply a bizarre scaling factor.
- The names of the bands are SR_B4, SR_B3, SR_B2, for the visible bands.
- For surface reflectance, we don’t need to worry about clouds.
- Because of the algorithms this data is created with, we don’t need to worry about creating a mean reflectance value. I won’t go into the details, but there are several reasons we don’t apply a mean.
Let’s look at how we apply this to the burn area over the Badger Fire from Exercise 2.1, question 1:
Example 3.5
- var dataset = ee.ImageCollection(‘LANDSAT/LC08/C02/T1_L2’)
- .filterDate(‘2023-07-01’, ‘2023-07-31’);
- // Applies scaling factors.
- function applyScaleFactors(image) {
- var opticalBands = image.select(‘SR_B.’).multiply(0.0000275).add(-0.2);
- return image.addBands(opticalBands, null, true);
- }
- dataset = dataset.map(applyScaleFactors);
- var visualization = {
- bands: [‘SR_B4’, ‘SR_B3’, ‘SR_B2’],
- min: 0.0,
- max: 0.3,
- };
- Map.setCenter(-114.1469, 42.2471, 10);
- Map.addLayer(dataset, visualization, ‘Landsat8 True Color’);
I’ve taken the scaling factor function directly from the example on the Earth Engine Catalog. You’re probably asking, why didn’t you just start with Landsat 8?!? It looks so much easier! Yes, it is. But look how much you grew by using the Sentinel imagery!