How to find a specific satellite image in Google Earth Engine

Let’s say there is a specific spatial event that you are interested in, for example the volcanic eruption in Tonga on January 15th in 2022.

For this example I will use data from publicly available polar orbiting satellites. These kind of satellites only take a single image of an area, and returns to that area perhaps a week later. So it is pretty unlikely that one of these satellites passed by at the exact time of the eruption, but I hope to find some before-and-after shots.

Now, you might have seen incredible images of events like the Tonga eruption. Here is a gif of the Japanese meteorological satellite Himawari 8:

From the Japanese National Institute of Informatics (link)

Meteorological satellites are geostationary, which means that they follow the rotation of the earth and therefore appear as if they are in a fixed position in the sky. Geostationary satellites are typically placed at an altitude of 35,800 kilometers, compared to only 500-600 kilometers of polar orbiting satellites. A satellite like the Himawari 8 can obtain images e.g. every 10 minutes, which is not what we will be dealing with now (sorry to disappoint).

The first thing to do in Google Earth Engine is to get the location of the event. In this case, it is the Hunga Tonga–Hunga Ha’apai volcano in the southern Pacific Ocean.

I create a geometry point at approximately (-175.39, -20.55), which roughly corresponds to where the volcano is located. I also create a polygon, which I will use to search for satellite images that overlap the area.

Point and polygon are in place

Now let’s hunt for satellite images!

The two Earth Observation satellite programmes that I am most familiar with are Landsat and Sentinel, so let’s focus on those.

I want to start by getting an idea of what the volcano looked like before the eruption. Let’s get a nice, cloud-free image from both satellite programmes.

Since Landsat 9 was launched in September 2021, I will use data from Landsat 8 instead for this pre-eruption image. I will set the cloud cover to a maximum of 20% to keep the image as clear as possible:

var l8_2021 = ee.ImageCollection('LANDSAT/LC08/C02/T1_TOA')
                  .filterDate('2021-01-01', '2022-01-01')
                  .filterBounds(geometry)
                  .filter(ee.Filter.lt('CLOUD_COVER_LAND', 20));
var l8_2021 = l8_2021.select(['B2','B3','B4'], ['B2','B3','B4'])
print(l8_2021, 'Landsat 8 (2021 average)');

var visualisation_l8 = {
  bands: ['B4', 'B3', 'B2'],
  min: 0.01,
  max: 0.2,
};

Map.addLayer(l8_2021.mean(), visualisation_l8, 'Landsat 8 2021')

Let’s do the same for Sentinel 2:

var s2_2021 = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2021-01-01', '2022-01-01')
                  .filterBounds(geometry)
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 5));

var s2_2021 = s2_2021.select(['B2','B3','B4'], ['B2','B3','B4'])
print(s2_2021, 'Sentinel 2 (2021 average)');

var visualisation_s2 = {
  min: 600,
  max: 2000,
  bands: ['B4', 'B3', 'B2'],
};

Map.addLayer(s2_2021.mean(), visualisation_s2, 'Sentinel 2 2021');

Note that I set the cloud cover to only 5% for the Sentinel 2 images. I allow myself to set it that low, because there are much more images available than in the Landsat collection.

Here are the resulting images:

The Landsat 8 (left) is comprised of the average of 3 images,
whereas the Sentinel 2 (right) is based on the average of 20 images.

From this comparison I think it is quite clear that the Sentinel 2 images have a higher spatial resolution than the Landsat 9 images.

Now we know the situation before the eruption. Next we need to find out if there was a satellite passing over the volcano on the day of the eruption, or at least closely after. First, let’s have a look at the newly launched Landsat 9:

var l9 = ee.ImageCollection('LANDSAT/LC09/C02/T1_TOA')
                  .filterDate('2022-01-15', '2022-01-20')
                  .filterBounds(geometry);
var l9 = l9.select(['B2','B3','B4'], ['B2','B3','B4'])
print(l9, 'Landsat 9 (eruption)');

From the print-statement I can see that there is a single image from Landsat 9 available, from January 17th:

Let’s do the same for Sentinel 2:

var s2 = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2022-01-15', '2022-01-20')
                  .filterBounds(geometry);

var s2 = s2.select(['B2','B3','B4'], ['B2','B3','B4'])
print(s2, 'Sentinel 2 (eruption)');

Here I also only get a single image, also from January 17th:

Luckily, both images are fairly cloud free, so we are able to see what remains of the island after the volcanic eruption two days prior:

Landsat 9 (left) and Sentinel 2 (right) on January 17 2022, two days after the eruption.

When comparing one of the images from January 17th to the average from 2021, it becomes evident that almost the entire island has disappeared as a result of the eruption:

I drew a polygon to make the difference more clear in the post-eruption image:

According to an estimate by NASA’s Earth Observatory, the eruption of Hunga Tonga–Hunga Ha’apai “released hundreds of times the equivalent mechanical energy of the Hiroshima nuclear explosion”. The plume from the eruption rose to 58 kilometers, making it the largest known volcanic plume measured by satellites.

You can find the Google Earth Engine code here.

SMAP soil moisture data now available at 10km resolution [Google Earth Engine]

I realise that I am a bit late to the party, but I just noticed the paper by Mladenova et al (2020) which outlines a new SMAP-based soil moisture dataset, where the resolution is increased from 0.25°x0.25° (roughly 27-28km at the equator) to 10km by “assimilating soil moisture retrievals from the Soil Moisture Active Passive (SMAP) mission into the USDA-FAS Palmer model“.

The researchers tested the new dataset in three areas: California (USA), the Western Cape (South Africa) and New South Wales (Australia), as all three areas experienced a drought since the SMAP satellite has been operational.

Study areas for the Mladenova et al (2020) paper.
Study area for the Mladenova et al (2020) paper

Soil moisture is an often overlooked dataset, especially due to the often terribly low spatial resolution. Perhaps this enhanced dataset will result in an increased usage of the data. The increased spatial resolution is, in my opinion, a big improvement:

0.25°x0.25° SMAP image compared to the 10x10km enhanced SMAP image.
Here I used Lithuania used an example, with data from 2017-2018.

As you can guess from the two images above, the enhanced dataset is available on Google Earth Engine.

If you want to see the (to my knowledge) only high-resolution soil moisture product in existence, check out the Dutch remote sensing company VanderSat. Their product is based on passive microwaves, and has a spatial resolution of just 100m.

The NASA-USDA soil moisture product generated from the SMAP-enhanced Palmer model can be found here: https://gimms.gsfc.nasa.gov/SMOS/SMAP/.

Time series for SSM, ET and precipitation for Myanmar [Google Earth Engine]

Back in December 2018, I was in Yangon (Myanmar) to help facilitate a workshop on satellite-based remote sensing to support water professionals in the country. For this workshop, I created a script to calculate a time series for surface soil moisture, evapotranspiration and precipitation.

The data for this script are:

The first thing to do is to load the shapefile of Myanmar, using the LSIB dataset. We need to deal with the auto generated system:index, which otherwise causes some visualisation issues later on. The auto generated ID for the Feature could be something like 000000000000000092.

My solution is a bit of a ‘hack’, because in this specific case, I only have a single Feature in my FeatureCollection, i.e. the shapefile of Myanmar. That means that I can use .first() and .set() to change the system:index of the single Feature. In this instance, I change the system:index to Myanmar.

// Select country from list
var countries = ee.FeatureCollection("USDOS/LSIB/2017");
var fc = countries.filter(ee.Filter.inList('COUNTRY_NA',['Burma']));
print('Old system:index', fc);

// Change system:index from the auto generated name
var feature = ee.Feature(fc.first().set('system:index', 'Myanmar')); // .first() is used as there is just one Feature
var Myanmar = ee.FeatureCollection(feature);
print('New system:index', Myanmar);

// Add Myanmar and center map
Map.centerObject(Myanmar);
Map.addLayer(Myanmar, {}, 'Myanmar shapefile');

This gives the following output:

Now I load the imageCollections of the three datasets.

// Filter datasets for bands and period
var SSM = ee.ImageCollection('NASA_USDA/HSL/SMAP_soil_moisture')
      .select('ssm')
      .filterBounds(Myanmar)
      .filterDate(startdate, enddate);

var ET = ee.ImageCollection('MODIS/006/MOD16A2')
      .select('ET')
      .filterBounds(Myanmar)
      .filterDate(startdate, enddate);
      
var PRCP = ee.ImageCollection('UCSB-CHG/CHIRPS/DAILY')
      .select('precipitation')
      .filterBounds(Myanmar)
      .filterDate(startdate, enddate);

Next, I create a single image consisting of the mean of each imageCollection for the given time period. This is purely for visualisation purposes as you will see next.

// Calculate means for visualisation
var SSM_mean = SSM.mean()
  .clip(Myanmar);

var ET_mean = ET.mean()
  .clip(Myanmar);

var PRCP_mean = PRCP.mean()
  .clip(Myanmar);

I set the visualisation parameters, and add the data to my map.

// Set visualisation parameters    
var soilMoistureVis = {
  min: 10.0,
  max: 28.0,
  palette: ['0300ff', '418504', 'efff07', 'efff07', 'ff0303'],
};

var evapotranspirationVis = {
  min: 0.0,
  max: 300.0,
  palette: [
    'ffffff', 'fcd163', '99b718', '66a000', '3e8601', '207401', '056201',
    '004c00', '011301'
  ],
};

var precipitationVis = {
  min: 1.0,
  max: 30.0,
  palette: ['#ffffcc','#a1dab4','#41b6c4','#2c7fb8','#253494'],
};

// Visualize maps
Map.addLayer(SSM_mean, soilMoistureVis, 'SSM');
Map.addLayer(ET_mean, evapotranspirationVis, 'ET');
Map.addLayer(PRCP_mean, precipitationVis, 'PRCP');

That gives us the following output:


Now what we really want is a graph showing the development of the data through time. Due to the limited memory capacity of Google Earth Engine, I can’t get more than 1.5 month of data, but that is enough to have a look at the onset of the monsoon season in May/June.

I use the ui.Chart.image.seriesByRegion where the arguments are: imageCollection, regions, reducer, band, scale, xProperty, seriesProperty (last four are optional). In addition, I use .setOptions() to add a main title and axis titles.

// Define graphs
var SSMGraph = ui.Chart.image.seriesByRegion({
  imageCollection: SSM, 
  regions: Myanmar, 
  band: 'ssm',
  reducer: ee.Reducer.mean(),
  scale: 25000, // 25km resolution
}).setOptions({
        title: 'Average soil moisture',
        hAxis: {title: 'Date'},
        vAxis: {title: 'Soil moisture [mm]'}});
        
var ETGraph = ui.Chart.image.seriesByRegion({
  imageCollection: ET, 
  regions: Myanmar, 
  band: 'ET',
  reducer: ee.Reducer.mean(),
  scale: 500, // 500m resolution
}).setOptions({
        title: 'Average evapotranspiration',
        hAxis: {title: 'Date'},
        vAxis: {title: 'Evapotranspiration [kg/m^2]'}});

var PRCPGraph = ui.Chart.image.seriesByRegion({
  imageCollection: PRCP, 
  regions: Myanmar, 
  band: 'precipitation',
  reducer: ee.Reducer.mean(),
  scale: 5000, // ~5km resolution
}).setOptions({
        title: 'Average precipitation',
        hAxis: {title: 'Date'},
        vAxis: {title: 'Precipitation [mm/day]'}});

// Display graphs
print(SSMGraph);
print(ETGraph);
print(PRCPGraph);

Which gives us the following output:

Note that “Myanmar” appears in the legend of the graph. If we had not changed the system:index in the beginning, you would instead see something like “000000000000000092“.

Link to the entire script: https://code.earthengine.google.com/e9042e4ac5fc902d1ff58da10856965f

Map examples

Here I will post examples of maps and other images that I have created. Feel free to ask about how certain elements were created, and I can make a post about it.

I will update this post when I have new content to show.

English and Danish version of a map showing the year of construction of buildings in the Danish capital, Copenhagen. Created using QGIS 3.10. Read about how I made the map here.
Geological formation that appeared in the VanderSat soil moisture data. Created using QGIS.

Calculate the percentage (%) cropland for each polygon in a shapefile [Google Earth Engine]

In a previous post, I use a similar script to calculate NDVI per polygon. However, this script includes a small difference which I have not yet myself understood why it was necessary.


As before, my shapefile consists of the administrative borders of Tanzania. I got this shapefile from the GADM database (version 3.6).

// Import administrative borders and select the gadm level 
// (0=country, 1=region, 2=district, 3=ward)
var gadm_level = '2';
var fc_oldID = ee.FeatureCollection('users/username/gadm36_TZA_'+gadm_level);
Map.addLayer(fc, {}, 'gadm TZA');
Map.centerObject(fc, 6);

Now I change the system:index of the FeatureCollection.

Why do this at all? I was struggling to produce a mean value of the COPERNICUS landcover data, after having created a successful script for the MODIS NDVI dataset. When I ran the script, no mean column was created. However, I was able to calculate mean values from a couple of polygons that I created myself within Google Earth Engine, so I knew the problem had to be in the shapefile somehow. So I started to investigate the gadm asset (i.e. the shapefile) and saw that Google Earth Engine automatically assigns some odd index names to each feature (e.g. 000000000001 or 0000000000a). I used the following steps described by the user Kuik on StackExchange to change the system:index:

// Change the system:index ("ID")

var idList = ee.List.sequence(0,fc_oldID.size().subtract(1));

var list = fc_oldID.toList(fc_oldID.size());

var fc = ee.FeatureCollection(idList.map(function(newSysIndex){
  var feat = ee.Feature(list.get(newSysIndex));
  var indexString = ee.Number(newSysIndex).format('%01d');
  return feat.set('system:index', indexString, 'ID', indexString);
}));

print('new system:index', fc);

Now you can see that the system:index has been changed:

Now I load in the landcover data. Here I use the 100m resolution, yearly Copernicus Global Land Cover dataset. Here we can find an image per year from 2015-2019. This dataset has a long list of interesting bands, e.g. a land cover classification with 23 categories, percentage of urban areas and percentage of cropland. In this script, I will look at the percentage of cropland (crops-coverfraction) for different administrative regions in Tanzania. Read more about the COPERNICUS dataset here.

// Load data

var year = '2015'
var lc = ee.Image("COPERNICUS/Landcover/100m/Proba-V-C3/Global/"+year)
      .select('crops-coverfraction')
      .clip(fc)

print('Image from '+year, lc)

print('Projection, crs, and crs_transform:', lc.projection());
print('Scale in meters:', lc.projection().nominalScale());

The last two lines is a good way to get a bit more information about your data if you, like I was, are looking for a reason why your script is misbehaving. It will give you the following output:

Now we can use the reducer reduceRegions. Again, remember that the scale should reflect the spatial resolution of the data. In the image above, we can see that the scale is 100 m.

// Reduce regions
var MeansOfFeatures= lc.reduceRegions({
  reducer: ee.Reducer.mean(),
  collection: fc,
  scale: 100,
})

print('gadm mean', ee.Feature(mean.first())) // Test if mean is produced

// Set visualization parameters    
var lcVis = {
  min: 0.0,
  max: 100.0,
  palette: [
    'ffffff', 'fcd163', '99b718', '66a000', '3e8601', '207401', '056201',
    '004c00', '011301'
  ],
};

Map.addLayer(lc, lcVis, 'crop cover (%) '+year)

// Export the featureCollection (table) as a CSV file
Export.table.toDrive({
  collection: MeansOfFeatures,
  description: 'proba_V_C3_TZA_adm'+gadm_level+'_'+year,
  folder: 'GEE output',
  fileFormat: 'CSV'
});

This should give you the following output:


You can find the full script here.

Calculate a mean NDVI value for each polygon in a shapefile [Google Earth Engine]

To run this script, you must have uploaded a shapefile to your Google Earth Engine account (an “Asset”). In my case, I am using the administrative borders for Tanzania. This dataset comes in three different levels. You can find that kind of data in the GADM database (version 3.6).

// Import administrative borders and select the gadm level 
// (0=country, 1=region, 2=district, 3=ward)
var gadm_level = '1';
var area = ee.FeatureCollection('users/username/gadm36_TZA_'+gadm_level);
print('area', area);
Map.centerObject(area, 6);
Map.addLayer(area, {}, 'area');

Next, you set the date of the period that you are interested in.

// Set the date
var year = '2019'
var startdate = year+'-01-01';
var enddate = year+'-12-31';

print(startdate)

Now I load in the NDVI data. Here I use the 250m resolution, 16-day Vegetation Indices product from MODIS (MOD13Q1.006). This dataset has a pre-calculated NDVI band, so I don’t have to first do the calculation to NDVI from the NIR/Red bands. However, in Earth Engine, the scale is 0.0001, so we need to correct that scale to get values in the order that we are used to deal with (from -1 to 1).

// MODIS NDVI scale is 0.0001, so we need to correct the scale
var NDVI = function(image) {
  return image.expression('float(b("NDVI")/10000)')
};

var modis = ee.ImageCollection('MODIS/006/MOD13Q1')
  .filterBounds(area)
  .filterDate(startdate,enddate);

// Now we select only the NDVI from the list of bands
var modisNDVI = modis.map(NDVI)
print('MODIS after .map', modisNDVI) 

Now we have the variable modisNDVI which is an ImageCollection. The ImageCollection contains a list of images that are within the start date and end date that we have specified in the beginning.

Now we’ll create a mosaic from this ImageCollection. You can read more about what a mosaic is here.

var mosaic = modisNDVI.mosaic()
  .clip(area);
print('mosaic', mosaic);

Another thing we could have done is to take the mean:

var mean = modisNDVI.mean()
  .clip(area);
print('mean', mean);

However, at least in my case, this takes much longer to compute.

Next, we set the visualisation parameters and add the mosaic to the map.

// Set visualisation parameters
var colour = {
 min: 0.0,
 max: 1.0,
 palette: [
 'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
 '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
 '012E01', '011D01', '011301'
 ],
};

Map.addLayer(mosaic, colour, 'spatial mosaic');

Now I use reduceRegions to get a value for each polygon (or region) in the shapefile. Remember that the scale should represent the spatial resolution of the dataset. Print out the result of the first polygon to make sure that you now have a column called “mean”.

// Add reducer output to the Features in the collection.
var MeansOfFeatures = mosaic.reduceRegions({
  collection: area,
  reducer: ee.Reducer.mean(),
  scale: 250,
});

print(ee.Feature(MeansOfFeatures.first()))

That should look something like this:

If you export the Feature now, you will notice that you get a column called “.geo“, which can become quite annoying if you are working with Excel, or just want your output to look nice. The .geo column contains the geometry, i.e. the XY coordinates, of each and every point of each and every polygon. This takes up a lot of space and is usually not data that you are interested in.
The user Rodrigo E. Principe gave a simple solution to get rid of the .geo-column in this post on StackExchange:

// Remove the content of .geo column (the geometry)
var featureCollection = MeansOfFeatures.map(function(feat){
  var nullfeat = ee.Feature(null)
  return nullfeat.copyProperties(feat)
})

Now we can export our output:

// Export the featureCollection (table) as a CSV file to your Google Drive account
Export.table.toDrive({
  collection: featureCollection,
  description: 'modis250_ndvi_TZA_adm'+gadm_level+'_'+year,
  folder: 'GEE output',
  fileFormat: 'CSV'
});

Your output will look something like this:

You can see the code in its entirety here.