1

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

8
  • Please, define what do you understand by a "pixel-wise comparison of the correlation between the datasets". Pixel to pixel comparison, for a same date, it has no sense because you will have one value for each dataset. So, Pearson's Correlation will reduce all the pixel values into a single statistic when you use your region (with almost 20,000,000 of pixels) for a same date.
    – xunilk
    Jul 23 at 14:06
  • Okay. So I basically am trying to find the degree of correlation between LST taken from MODIS sensor and Landsat LST sensor. As you said, using reduceRegion() will reduce all the pixels into a single statistic, but I am unsure how that actually works with regard to the Pearson's Correlation. I understand I have to create an Image collection such that each Image has two bands for which I want the correlation coefficient. Will the output be Images with each pixel containing correlation values? If yes will I have to plot that with a reducer.mean() ? Thanks Jul 24 at 12:25
  • If you have any other suggestions as to how I can visualize the degree of correlation between the LST values taken by MODIS sensors and Landsat over the roi ? (If not the Pearson's correlation) Jul 24 at 12:47
  • First of all, why do you want to obtain the correlation between Modis and Landsat? For the daily basis of Modis? If you work on a large scale, it probably works; even though the images have masked regions (the series will be unbalanced with respect to Landsat). At small scales it is shortsighted to detect changes due to its low resolution.
    – xunilk
    Jul 24 at 14:39
  • Well, the region I am working with is quite remote and does not have many weather stations. There has been a lot of research done comparing satellite data with weather station data in terms of LST, however, there are not many researches comparing the correlation or congruence between different satellite sensors. I am trying to see whether an upscaled sample of MODIS will can be used as a "high temporal resolution" alternative to Landsat. (Even though Landsat resolution is better, the MODIS has the advantage of daily coverage). Jul 24 at 15:26

1 Answer 1

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.

enter image description here

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.

enter image description here

10
  • Thankyou very much for your answer. I guess I could also chart this list for better visualization. And Instead of Pearson's correlation I could just do linear regression reducer. Jul 25 at 12:19
  • You're welcome. Yes, you can. However, when you map for different dates remember that Landsat series have images with some repeated dates and incomplete coverage (Modis is on a global scale). For this reason both series are unbalanced.
    – xunilk
    Jul 25 at 13:28
  • Hi @xunilk. I know this is a new question, and I have also posted it separately, but I would appreciate it if you could help me out in this. So I plotted the mean_LST layer and Landsat LST values on a scatter plot, and even after bitmasking, I am getting many anomalous values in the mean_LST layer such as temperature values of -143 C. What could be the reason for these outliers and how can I remove if possible ? Thanks in advance. Here is the link for the code: code.earthengine.google.com/faff50d540af6a40386d9a4d80d562a5 If it takes too long to process just reduce the time period. Aug 2 at 0:17
  • I will take a look. First of all, you need to corroborate that each paired value correspond effectively to the same date and point sampled. I ran the script and it seems that it is not true.
    – xunilk
    Aug 2 at 2:09
  • How do you mean ? As far as I can see the size of the merged collection looks okay. The Landsat image collection has 604 elements, and the merged collection has 514 elements. Instead of sampling points, I am using a reducer to calculate the total mean of each image and then plotting on the scatter plot. Aug 2 at 20:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.