# Temporal Gap-Filling with Linear Interpolation in GEE

Many applications require replacing missing pixels in an image with an interpolated value from its temporal neighbours. This gap-filling technique is used in several applications, including:

• Replacing Cloudy Pixels: You may want to fill gap in an image with the best-estimated value from the before and after cloud-free pixel.
• Estimating Intermediate Values: You can use this technique to compute an image for a previously unknown time-step. If you had population rasters at 2 different years and want to compute a population raster for an intermediate year using pixel-wise linear interpolation.
• Preparing Data for Regression: All of your independent variables may not be available at the same temporal resolution. You can harmonize various dataset by generating interpolated raster at uniform or fixed time-steps.

Google Earth Engine can be used effectively for gap-filling time-series datasets. While the logic for linear interpolation is fairly straightforward, data preparation for this in GEE can be quite challenging. It involves use of Joins, Masks and Advanced Filters. This post explains the steps with code snippets and builds a fully functioning script that can be applied on any time-series data.

You can also watch my video at Geo for Good 2022 Conference where I covered the same topic.

### Overview Of Steps

Interpolating an Image Collection in GEE involves the following steps

• Adding a time-band to each image
• Joining the image collection with itself to find images before and after each image
• Processing the result to create interpolated images and replacing masked pixels with interpolated values.
• Visualizing/Exporting the results

### Example 1: Gap-Filling MAsked NDVI Pixels

We will see how to apply this technique on a Sentinel-2 Time-Series and replace cloudy pixels with interpolated values from images taken before and after it. The image below shows the result of interpolation.

We start by preparing a time-series for 1-year by filtering the Sentinel-2 Level-2A collection and applying a cloud-masking function. The code below uses a polygon for a farm as a geometry. You may replace it with the geometry for your area of interest.

```var geometry = ee.Geometry.Polygon([[
[82.60642647743225, 27.16350437805251],
[82.60984897613525, 27.1618529901377],
[82.61088967323303, 27.163695288375266],
[82.60757446289062, 27.16517483230927]
]]);
Map.centerObject(geometry)
var s2 = ee.ImageCollection("COPERNICUS/S2_SR");
var filtered = s2
.filter(ee.Filter.date('2019-01-01', '2020-01-01'))
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
.filter(ee.Filter.bounds(geometry))
// Write a function for Cloud masking
var cloudProb = image.select('MSK_CLDPRB');
var snowProb = image.select('MSK_SNWPRB');
var cloud = cloudProb.lt(5);
var snow = snowProb.lt(5);
var scl = image.select('SCL');
var cirrus = scl.eq(10); // 10 = cirrus
// Cloud probability less than 5% or cloud shadow classification
.select("B.*")
.copyProperties(image, ["system:time_start"]);
}

```

Add a band containing timestamp to each image. This will be used to do pixel-wise interpolation later. The key here is to mask the timestamp band with the same mask as the image. If your image has more than 1-band, the code below will use the mask from the first band.

```var filtered = filtered.map(function(image) {
})
```

Next we specify the time-window. This will determine how much backward and forward are we willing to look for an unmasked pixel in the time-series. The filters in subsequent steps need the time difference in milliseconds, so we convert days to milliseconds.

```var days = 30
var millis = ee.Number(days).multiply(1000*60*60*24)
```

For each image in the collection, we need to find all images before and after the specified time-window. This is accomplished using Joins. We need to do 2 joins

• Join 1: Join the collection with itself to find all images before each image
• Join 2: Join the collection with itself to find all images after each image

We first define the filters needed for the join. Define a maxDifference filter to find all images within the specified days.

```var maxDiffFilter = ee.Filter.maxDifference({
difference: millis,
leftField: 'system:time_start',
rightField: 'system:time_start'
})
```

We need 2 more filters.

• A `lessThanOrEquals` filter to find all images after a given image. This will compare the given image’s timestamp against other images’ timestamps
• A `greaterThanOrEquals` filter to find all images before a given image. This will compare the given image’s timestamp against other images’ timestamps
```var lessEqFilter = ee.Filter.lessThanOrEquals({
leftField: 'system:time_start',
rightField: 'system:time_start'
})
var greaterEqFilter = ee.Filter.greaterThanOrEquals({
leftField: 'system:time_start',
rightField: 'system:time_start'
})
```

Now we apply the joins. We join the image collection with itself and match before/after images for each image in the collection that meet our criteria.

Join1 is to find all images after each image that are winthin the time-window. We have 2 filters already defined for these two conditions. We combine them and apply the join. This join will find all images after, sorted in descending order. This will gives us resulting images such that closest image is last.

```var filter1 = ee.Filter.and(maxDiffFilter, lessEqFilter)

var join1 = ee.Join.saveAll({
matchesKey: 'after',
ordering: 'system:time_start',
ascending: false})

var join1Result = join1.apply({
primary: filtered,
secondary: filtered,
condition: filter1
})
```

Do the Join2 now to match all images within the time-window that come before each image. This join will find all images before, sorted in ascending order. This will gives us images so that closest is last.

```var filter2 = ee.Filter.and(maxDiffFilter, greaterEqFilter)

var join2 = ee.Join.saveAll({
matchesKey: 'before',
ordering: 'system:time_start',
ascending: true})

var join2Result = join2.apply({
primary: join1Result,
secondary: join1Result,
condition: filter2
})
```

The collection is now ready for interpolation. We write a function that will be used to interpolate all images. This function takes an image and replaces the masked pixels with the interpolated value from before and after images.

```Interpolation Formula
y = y1 + (y2-y1)*((t – t1) / (t2 – t1))
y = interpolated image
y1 = before image
y2 = after image
t = interpolation timestamp
t1 = before image timestamp
t2 = after image timestamp
```

Tip:

If you wanted a simple average of before and after images, you can set the value of `timeRatio` to 0.5.

```var interpolateImages = function(image) {
var image = ee.Image(image)
// We get the list of before and after images from the image property
// Mosaic the images so we a before and after image with the closest unmasked pixel
var beforeImages = ee.List(image.get('before'))
var beforeMosaic = ee.ImageCollection.fromImages(beforeImages).mosaic()
var afterImages = ee.List(image.get('after'))
var afterMosaic = ee.ImageCollection.fromImages(afterImages).mosaic()
// Get image with before and after times
var t1 = beforeMosaic.select('timestamp').rename('t1')
var t2 = afterMosaic.select('timestamp').rename('t2')
var timeImage = ee.Image.cat([t1, t2, t])
var timeRatio = timeImage.expression('(t - t1) / (t2 - t1)', {
't': timeImage.select('t'),
't1': timeImage.select('t1'),
't2': timeImage.select('t2'),
})
// Compute an image with the interpolated image y
var interpolated = beforeMosaic
// Replace the masked pixels in the current image with the average value
return result.copyProperties(image, ['system:time_start'])
}
```

Now `map()` the function to interpolate all images in the collection. The `interpolatedCol` is a continuous time-series with cloudy pixels replaced by interpolated values from closest cloud-free pixels.

```var interpolatedCol = ee.ImageCollection(
join2Result.map(interpolateImages))
```

Here’s the full script along with code for computing NDVI and creating an animation showing the gap-filled collection. https://code.earthengine.google.co.in/536d8fc3b031365e6d4f75d38131682e

### Example 2: Interpolating Population Rasters

Another useful application of this technique is to generate data at missing time intervals. We will take the GHSL: Global Human Settlement Layers Population Grid at year 2000 and 2015 and interpolate it to create population grids from year 2005 and 2010. Below image shows the results we are trying to achieve.

We can follow the same process as the previous example with a small change. The previous example already had images at all the time-steps that we needed and we only had to fill the masked pixels with interpolated values. Here, we need prepare our collection so it has empty (fully masked) images for the years that we want to generate new population grids. The code will then fill all the pixels with interpolated values.

```var pop2000 = ee.Image('JRC/GHSL/P2016/POP_GPW_GLOBE_V1/2000')
var pop2015 = ee.Image('JRC/GHSL/P2016/POP_GPW_GLOBE_V1/2015')
var ghsl = ee.ImageCollection.fromImages([pop2000, pop2015])
var yearsToInterpolate = ee.List([2005, 2010]);
var initImages = yearsToInterpolate.map(function(year) {
var image = ee.Image().rename('population_count').toFloat().set({
'system:index': ee.Number(year).format(),
'system:time_start': ee.Date.fromYMD(year, 1, 1). millis(),
})
return image
})
var initCol = ee.ImageCollection.fromImages(initImages)
var mergedCol = ghsl.merge(initCol)
var mergedCol = mergedCol.map(function(image) {
})
```

Next, we need to choose a time-window. We have population images 15 years apart, so we set time-window to be `15 x 365` days.

```var days = 15*365
var millis = ee.Number(days).multiply(1000*60*60*24)
```

Rest of the code for interpolation remains the same. We just need to change the variable name for the collection in join1.

```var maxDiffFilter = ee.Filter.maxDifference({
difference: millis,
leftField: 'system:time_start',
rightField: 'system:time_start'
})
var lessEqFilter = ee.Filter.lessThanOrEquals({
leftField: 'system:time_start',
rightField: 'system:time_start'
})
var greaterEqFilter = ee.Filter.greaterThanOrEquals({
leftField: 'system:time_start',
rightField: 'system:time_start'
})
var filter1 = ee.Filter.and(maxDiffFilter, lessEqFilter)
var join1 = ee.Join.saveAll({
matchesKey: 'after',
ordering: 'system:time_start',
ascending: false})

var join1Result = join1.apply({
primary: mergedCol,
secondary: mergedCol,
condition: filter1
})
var filter2 = ee.Filter.and(maxDiffFilter, greaterEqFilter)
var join2 = ee.Join.saveAll({
matchesKey: 'before',
ordering: 'system:time_start',
ascending: true})

var join2Result = join2.apply({
primary: join1Result,
secondary: join1Result,
condition: filter2
})
var interpolateImages = function(image) {
var image = ee.Image(image)
var beforeImages = ee.List(image.get('before'))
var beforeMosaic = ee.ImageCollection.fromImages(beforeImages).mosaic()
var afterImages = ee.List(image.get('after'))
var afterMosaic = ee.ImageCollection.fromImages(afterImages).mosaic()
var t1 = beforeMosaic.select('timestamp').rename('t1')
var t2 = afterMosaic.select('timestamp').rename('t2')
var timeImage = ee.Image.cat([t1, t2, t])
var timeRatio = timeImage.expression('(t - t1) / (t2 - t1)', {
't': timeImage.select('t'),
't1': timeImage.select('t1'),
't2': timeImage.select('t2'),
})
var interpolated = beforeMosaic
return result.copyProperties(image, ['system:time_start'])
}
var interpolatedCol = ee.ImageCollection(join2Result.map(interpolateImages))

```

The interpolated collection is ready. We can visualize the images by adding each of them to the Map. Here we are using the boundary for City of Johannesburg, South Africa and clipping each image to it.

```// Add a 'year' property to each image so we can filter them easily
var finalCol = interpolatedCol.map(function(image) {
var date = ee.Date(image.get('system:time_start'))
return image.set('year', date.get('year'))
})
// Clip and Visualize the Images for a City
var johannesburg = ee.FeatureCollection('users/ujavalgandhi/public/COJ_Boundary')
var geometry = johannesburg.geometry()
// Select the population band and clip to the geometry
var pop2000Image = pop2000.select('population_count').clip(geometry)
var pop2015Image = pop2015.select('population_count').clip(geometry)
var pop2005Image = ee.Image(finalCol.filter(ee.Filter.eq('year', 2005)).first())
.select('population_count').clip(geometry)
var pop2010Image = ee.Image(finalCol.filter(ee.Filter.eq('year', 2010)).first())
.select('population_count').clip(geometry)
// Visualize on the Map
var palette = ['#fef0d9','#fdcc8a','#fc8d59','#e34a33','#b30000']
var popVisParams = {min:30, max:1500, palette: palette}
Map.centerObject(geometry)
```

You will find the full script for this example at https://code.earthengine.google.co.in/612e5b4dbc08995c9455a5fc7e7b1f55

This was one of the most frequently asked questions by many users and there was no open-source implementation that was readily available to use. I spent a significant amount of time developing this code and writing this post to help researchers apply this to their work. Do leave a comment if you found it useful or have ideas to improve the code.

If you are new to Earth Engine and want to master it, check out my course End-to-End Google Earth Engine.

1. Cindy Viviers says:

Wow what a rare and valuable piece of code.

How would I adjust the code of example 2 to create empty (fully masked images) for specific dates where data is missing?
Something like var monthsToInterpolate = ee.List([20020601,20020701,20030401]);

And then interpolate using the average of the previous month and sequential month’s reading? Or something like that.

cant seem to get this right.
Any insight will be greatly appreciated.

• ujaval says:

Thanks Cidy. Glad you found it useful. Below is a snippet to create images for the specific dates.

var monthsToInterpolate = ee.List(
[ee.Date.fromYMD(2002, 6, 1),
ee.Date.fromYMD(2002, 7, 1),
ee.Date.fromYMD(2003, 4, 1)])

var initImages = monthsToInterpolate.map(function(date) {
var year = ee.Date(date).get(‘year’).format();
var month = ee.Date(date).get(‘month’).format(‘%02d’);
var imageId = year.cat(month);
var image = ee.Image().rename(‘monnthlyimage’).toFloat().set({
‘system:index’: imageId,
‘system:time_start’: ee.Date(date). millis(),
})
return image
})
var initCol = ee.ImageCollection.fromImages(initImages)
print(initCol)

2. Benazir Meerasha says:

How to apply the above coding for my missing MODIS time series NDVI data.

3. AFEWORK MEKEBERIAW says:

Big thanks and have been following you for the last two years and your contribution is huge. Godbless you and your family

4. Christhian Santana Cunha says:

Thanks for so much! Your code is very useful for us to apply in many regions of Brazil where the incidence of clouds is high.

5. Enes says:

Thanks! It is very useful code.
I have a question related to example1 regarding the first and last images which have some nan values due to the lack of information before or after. How can I handle this problem.

• ujaval says:

The only solution is to increase the time period so there is an unmasked pixel before/after the image within that period. But increasing it too much can lead to filling hte pixel values from a different season. So you’ll have to use an appropriate value for your use case.

• Enes says:

Thanks for your response, I appreciate that.
Also, I have another question. Is it normal to be more missing pixels after interpolation of the first image compared to the original image?

• ujaval says:

No. The first image should have the same or fewer missing pixels.

• Enes says:

Yes, thanks. I was doing a mistake. Now it is working,

6. Mauricio Jean Freisleben says:

Hello! thanks for the tutorial. amazing.
However, i’m having some kind of problem here. Basically when i select a random image from the filteredcollection it gives me images with null values (holes), where it should be filled, right?
thanks!!!

• ujaval says:

Likely there were no cloud-free pixels in the range you selected. If you selected 30-day time-period and there were no cloud-free pixels during that time, you will still get holes.