Pixel-Wise comparison of MODIS and Landsat LST Values in GEE

Question

I am trying to find the correlation between the LST values obtained from MODIS Terra and Aqua, and Landsat 8 over a common roi. I have performed the following steps until now:

1. I took two image collections (Terra Day and Aqua Day) and used an inner join to filter only the images taken on the same date,
2. I upscaled the image collection using bilinear interpolation to a resolution of 100 to compare it with the Landsat images.
3. I created a new image collection containing the average LST values of each day in the time period by using the formula (TerraD + AquaD / 2). This was performed on the upscaled images.

I have also imported and clipped the Landsat 8 LST dataset. I now want to calculate the correlation between the mean daily LST dataset created from MODIS, and the Landsat LST collection. However, to calculate the Pearson's Correlation, I will have to use reduce region which will reduce all the pixel values into a single statistic when I actually want a pixel-wise comparison of the correlation between the datasets. How may I proceed with this? I also understand that I will have to compare images taken on the same date so I will have to repeat the filtration process for the mean LST and Landsat Dataset.

My code is given below:

var terraD = ee.ImageCollection('MODIS/061/MOD11A1')
  .filterDate('2022-01-01', '2023-01-01').select('LST_Day_1km')
  .filterBounds(geometry)

var aquaD = ee.ImageCollection('MODIS/061/MYD11A1')
  .filterDate('2022-01-01', '2023-01-01')
  .select('LST_Day_1km')
  .filterBounds(geometry);
  
var landsatD = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
  .filterDate('2022-01-01', '2023-01-01')
  .select('ST_B10')
  .filterBounds(geometry);
  
var landSurfaceTemperatureVis = {
  min: 13000.0,
  max: 16500.0,
  palette: [
    '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003'
  ],
};

// Function to clip each image in the ImageCollection to the ROI
var clipToROI = function(image) {
  return image.clip(geometry);
};

var clipTerra = terraD.map(clipToROI)
Map.addLayer(clipTerra, landSurfaceTemperatureVis, 'TerraD')

var clipAqua = aquaD.map(clipToROI)
Map.addLayer(clipAqua, landSurfaceTemperatureVis, 'AquaD')

var clipLandsat = landsatD.map(clipToROI)
Map.addLayer(clipLandsat)

var terraDayCount = clipTerra.size().getInfo();
if (terraDayCount > 0) {
  print('MODIS Terra daytime data is available. Count:', terraDayCount);
} else {
  print('MODIS Terra daytime data is unavailable for the specified date range.');
}
//////////UPSCALE////////////////////

// Function to upscale an image using bilinear interpolation
var upscaleBilinear = function(image) {
  return image.resample('bilinear').reproject({
    crs: image.projection(),
    scale: 100  // Set the desired scale (resolution)
  });
};

// Apply bilinear interpolation to the Terra and Aqua datasets
var bilinearTerra = clipTerra.map(upscaleBilinear);
var bilinearAqua = clipAqua.map(upscaleBilinear);
print(bilinearTerra)

// Add the upscaled Terra and Aqua layers to the map with the specified visualization
Map.addLayer(bilinearTerra, landSurfaceTemperatureVis, 'MODIS Terra (Upscaled)');
Map.addLayer(bilinearAqua, landSurfaceTemperatureVis, 'MODIS Aqua (Upscaled)');

// Join Terra and Aqua images based on acquisition date
var join = ee.Join.inner().apply({
  primary: bilinearTerra,
  secondary: bilinearAqua,
  condition: ee.Filter.equals({
    leftField: 'system:time_start',
    rightField: 'system:time_start'
  })
});

//////////////////////MEAN////////////////////////

// Function to calculate the mean of Terra and Aqua images
var calculateMean = function(image) {
  // Get the Terra and Aqua images
  var terraImage = ee.Image(image.get('primary'));
  var aquaImage = ee.Image(image.get('secondary'));
  
  // Calculate the mean of Terra and Aqua images
  var meanImage = (terraImage.add(aquaImage)).divide(2).rename('mean_LST');
  
  // Return the mean image with the acquisition date
  return meanImage.set('system:time_start', terraImage.get('system:time_start'));
};

// Apply the calculateMean function to the joined ImageCollection
var meanCollection = ee.ImageCollection(join.map(calculateMean));

var first = meanCollection.first()
Map.addLayer(meanCollection,  landSurfaceTemperatureVis, 'mean' )

var matchedCount = meanCollection.size().getInfo();
if (matchedCount > 0) {
  print('Matching Terra and Aqua LST images found. Count:', matchedCount);
} else {
  print('No matching Terra and Aqua LST images found.');
}

print(meanCollection)
print(clipTerra)
print(clipAqua)
print(clipLandsat)

For viewing in GEE: Link

Answer 1

First of all, I could detect that "daily" Modis series have images with masked areas in considered dates range. So, as it can be observed at following picture, I only chose a little area (circled in red) of your region (with 9879 pixels) where it is was guaranteed not masked areas in Modis first image. I repeated the filtration process for both Modis and Landsat datasets and, 337 elements of initial meanCollection were converted in a Image Collection with same dates of Landsat collection (only 20 elements for each one). Datasets were both scaled as suggested by metadata information and temperatures are both converted in Celsius degrees.

Complete code looks as follows:

var pt = ee.Geometry.Point(75.858606, 32.407948);

Map.centerObject(pt, 16);

var terraD = ee.ImageCollection('MODIS/061/MOD11A1')
  .filterDate('2022-01-01', '2023-01-01').select('LST_Day_1km')
  .filterBounds(geometry);

var aquaD = ee.ImageCollection('MODIS/061/MYD11A1')
  .filterDate('2022-01-01', '2023-01-01')
  .select('LST_Day_1km')
  .filterBounds(geometry);
  
var landsatD = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
  .filterDate('2022-01-01', '2023-01-01')
  .select('ST_B10')
  .filterBounds(geometry)
  .map(function (img){
    return img.multiply(0.00341802).add(149).subtract(273.15)
              .set("system:time_start", img.get("system:time_start"));
  });

var landSurfaceTemperatureVis = {
  min: 13000.0,
  max: 16500.0,
  palette: [
    '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003'
  ],
};

// Function to clip each image in the ImageCollection to the ROI
var clipToROI = function(image) {
  return image.clip(geometry);
};

var clipTerra = terraD.map(clipToROI);
//Map.addLayer(clipTerra.first(), landSurfaceTemperatureVis, 'TerraD');

var clipAqua = aquaD.map(clipToROI);
//Map.addLayer(clipAqua.first(), landSurfaceTemperatureVis, 'AquaD');

var clipLandsat = landsatD.map(clipToROI);

//////////UPSCALE////////////////////

// Function to upscale an image using bilinear interpolation
var upscaleBilinear = function(image) {
  return image.resample('bilinear').reproject({
    crs: image.projection(),
    scale: 100  // Set the desired scale (resolution)
  });
};

// Apply bilinear interpolation to the Terra and Aqua datasets
var bilinearTerra = clipTerra.map(upscaleBilinear);
var bilinearAqua = clipAqua.map(upscaleBilinear);

//print("bilinearTerra", bilinearTerra);
//print("bilinearAqua", bilinearAqua);

// Add the upscaled Terra and Aqua layers to the map with the specified visualization
//Map.addLayer(bilinearTerra.first(), landSurfaceTemperatureVis, 'MODIS Terra (Upscaled)');
//Map.addLayer(bilinearAqua.first(), landSurfaceTemperatureVis, 'MODIS Aqua (Upscaled)');

// Join Terra and Aqua images based on acquisition date
var join = ee.Join.inner().apply({
  primary: bilinearTerra,
  secondary: bilinearAqua,
  condition: ee.Filter.equals({
    leftField: 'system:time_start',
    rightField: 'system:time_start'
  })
});

//////////////////////MEAN////////////////////////

// Function to calculate the mean of Terra and Aqua images
var calculateMean = function(image) {
  // Get the Terra and Aqua images
  var terraImage = ee.Image(image.get('primary'));
  var aquaImage = ee.Image(image.get('secondary'));
  
  // Calculate the mean of Terra and Aqua images
  var meanImage = terraImage.add(aquaImage)
    .divide(2)
    .multiply(0.02)
    .subtract(273.15)
    .rename('mean_LST');
  
  // Return the mean image with the acquisition date
  return meanImage.set('system:time_start', terraImage.get('system:time_start'));
};

// Apply the calculateMean function to the joined ImageCollection
var meanCollection = ee.ImageCollection(join.map(calculateMean));

print("meanCollection", meanCollection);

//print(meanCollection.first().projection().nominalScale());

var dates_meanCollection = meanCollection.aggregate_array("system:time_start").map(function (ele){
  
  return ee.Date(ele).format().slice(0,10);
  
}).distinct();

print("dates_meanCollection (distinct)", dates_meanCollection.sort());

var inner_landsat = clipLandsat.toList(clipLandsat.size()).map(function (ele){
  
  var date = ee.Date(ee.Image(ele).get('system:time_start')).format().slice(0,10);
  
  return ee.Algorithms.If(dates_meanCollection.contains(date), ele, 0);
  
}, true).removeAll([0]);

print("inner_landsat", inner_landsat);

var proj_modis = ee.ImageCollection(meanCollection)
  .first()
  .projection();

var new_modis_projection = ee.ImageCollection(meanCollection)
  .first()
  .reproject('EPSG:4326', null, proj_modis.nominalScale())
  .projection();

var grid_modis = geometry.coveringGrid(new_modis_projection);

print("grid_modis size", grid_modis.size());

Map.addLayer(clipLandsat.first(), {}, 'landsat8');

Map.addLayer(grid_modis, {}, 'grid_modis');

var dates_inner_landsat = ee.ImageCollection(inner_landsat).aggregate_array("system:time_start").map(function (ele){
  
  return ee.Date(ele).format().slice(0,10);
  
}).distinct();

var inner_meanCollection = meanCollection.toList(meanCollection.size()).map(function (ele){
  
  var date = ee.Date(ee.Image(ele).get('system:time_start')).format().slice(0,10);
  
  return ee.Algorithms.If(dates_inner_landsat.contains(date), ele, 0);
  
}, true).removeAll([0]);

print("inner_meanCollection", inner_meanCollection);


var temp_modis = grid_modis.toList(grid_modis.size()).map(function (ele) {
  
  var modis_temp = ee.ImageCollection(inner_meanCollection).first().sample({
    region: ee.Feature(ele).geometry(), 
    scale: 100, 
    projection: 'EPSG:4326'});

  return modis_temp.first().get('mean_LST');
  
});

print ("temp modis size", temp_modis.size());
print ("temp_modis date", ee.Date(ee.Image(inner_meanCollection.get(0)).get('system:time_start')).format().slice(0,10));

print("temp modis first image", ee.List(temp_modis).slice(0, 10));

var temp_landsat = grid_modis.toList(grid_modis.size()).map(function (ele) {
  
  var landsat_temp = ee.ImageCollection(clipLandsat).first().sample({
    region: ee.Feature(ele).geometry(), 
    scale: 100, 
    projection: 'EPSG:4326'});

  return landsat_temp.first().get('ST_B10');
  
});

print ("temp_landsat size", temp_landsat.size());
print ("temp_landsat date", ee.Date(ee.Image(inner_landsat.get(0)).get('system:time_start')).format().slice(0,10));

print("temp_landsat first image", ee.List(temp_landsat).slice(0, 10));

/////////////////Correlation/////////////////

After running it at GEE code editor, I got following result. It can be observed that Modis probably overestimates the temperatures, pixel to pixel, for the considered area. So, as I only took 10 pixels at the picture (for space reasons related to comparison), you can get an image for all area by calculating the differences pixel to pixel for corroborating that tendency. Afterward, map for other dates and if the tendency is corroborated it could be established a linear correlation for both datasets.