Null values after reduceRegion() and reducer.group

Question

I'm wondering why I get some null values after .reduceRegion() and reducer.group().

I'm extracting daily surface temperature values for many lakes (~400,000) from Landsat images during 1985-2021. The traditional way to do this is performing .reduceRegions() and reducer.mean() to all lake shapefiles. However, this method would be incredibly time-comsuming because of the complex boundary of lakes. For example, it took me 10~50 minutes to extract surface temperature for 130,000 lakes on one day. So I decided to first convert lake shapefiles to raster, with a new band called "Hylak_id" added to Landsat images, which indicates the lake id of each pixel. Then I performed zonal statistics according to this tutorial (https://developers.google.com/earth-engine/guides/reducers_grouping). But the result is null for many lakes. In the code below I presented two examples: tmpid = 1 (valid surface temperature) and tmpid = 1403425 (null). The weird thing is, if I computed .reduceRegions with the geometry of the lake 1403425, a valid surface temperature could be computed. I'm new to GEE and very confused about this.

Code editor link (I have shared my data): https://code.earthengine.google.com/dadb0c84bc013c8e8fd23c40f530d07a

// Two examples
var tmpid = 1403425
// var tmpid = 1

// Convert lake shapefile to raster
var lake_img = lakelist.filterMetadata('Hylak_id','equals',tmpid).reduceToImage({
  properties: ['Hylak_id'],
  reducer: ee.Reducer.first()
}).rename('Hylak_id')
var lake_mask = lake_img.gt(0)

// Read Landsat 5
var l5t1 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2').select(["QA_PIXEL", "ST_B6"])

// Create mask
var water = ee.Image('JRC/GSW1_4/GlobalSurfaceWater').select("occurrence").gte(90)
var water = water.updateMask(water.neq(0))

// A function to remove cloud, snow, and non-lake pixels, and add Hylak_id
function prep(img) {
  // remove cloud/snow and assign Hylak_id
  var systime = img.get('system:time_start')
  var qa = img.select(["QA_PIXEL"])
  var cloud1 = qa.bitwiseAnd(2).eq(0)  
  var cloud3 = qa.bitwiseAnd(8).eq(0) 
  var cloudshadow = qa.bitwiseAnd(16).eq(0) 
  var snow = qa.bitwiseAnd(32).eq(0)  
  var cloud_confid = qa.rightShift(8).bitwiseAnd(3).lt(2)  
  var cloudsh_confid = qa.rightShift(10).bitwiseAnd(3).lt(2) 
  var snow_confid = qa.rightShift(12).bitwiseAnd(3).lt(2)  
  var cirrus_confid = qa.rightShift(14).bitwiseAnd(3).lt(2)  
  var updated = (img.updateMask(cloud1)
  .updateMask(cloud3)
  .updateMask(snow)
  .updateMask(cloudshadow)
  .updateMask(cloud_confid)
  .updateMask(snow_confid)
  .updateMask(cirrus_confid)
  .updateMask(cloudsh_confid)
  .updateMask(water)
  .updateMask(lake_mask)
  ).addBands(lake_img)
  return updated
}

// Mosaic all images in the same day
function mosaicByDate(imcol){
  // imcol: An image collection
  // returns: An image collection
  var imlist = imcol.toList(imcol.size())
  var unique_dates = imlist.map(function(im){
    return ee.Image(im).date().format("YYYY-MM-dd")
  }).distinct()
  var mosaic_imlist = unique_dates.map(function(d){
    d = ee.Date(d)
    var im = imcol
      .filterDate(d, d.advance(1, "day"))
      .mosaic()  
    return im.set(
        "system:time_start", d.millis(), 
        "system:id", d.format("YYYY-MM-dd"))
  })
  return ee.ImageCollection(mosaic_imlist)
}

// Compute lake surface water temperature
function cal_lswt(img) {
  var lswt = img.select(['ST_B6']).multiply(0.00341802).add(149).add(-273.15).rename('surface_temp')
  return img.addBands(lswt)
}

// Choose one day for tests
var lswt = mosaicByDate(l5t1.filterDate('2002-06-10','2002-06-11').map(prep)).map(cal_lswt).select(['surface_temp', 'Hylak_id'])
var tmp_img = ee.Image(lswt.toList(lswt.size()).get(0))

var result = tmp_img.reduceRegion({
  reducer: ee.Reducer.mean().repeat(1).group({
    groupField: 1,
    groupName: 'Hylak_id'
  }),
  geometry: ee.Geometry.BBox(-180, -90, 180, 90), // I use a global geometry because I'd like to compute zonal statistics around the globe.
  crs: 'EPSG:4326',
  scale: 30,
  bestEffort: true,
  tileScale: 16,})

// Result using the reducer.group()
print(result)

// Result using the traditional way
print(tmp_img.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: lakelist.filterMetadata('Hylak_id','equals',tmpid),
  scale: 30,
  bestEffort: true,
  tileScale: 16
}))

// Show images
Map.addLayer(lake_img, null, 'Raster of lake polygon')
var vs_params =  {
  min: 2,
  max: 49,
  palette: [
    '000080', '0000d9', '4000ff', '8000ff', '0080ff', '00ffff', '00ff80',
    '80ff00', 'daff00', 'ffff00', 'fff500', 'ffda00', 'ffb000', 'ffa400',
    'ff4f00', 'ff2500', 'ff0a00', 'ff00ff'
  ]
}
Map.addLayer(tmp_img.select('surface_temp'), vs_params, 'Landsat 5')
var styling={color:'red', fillColor:'00000000'}
Map.addLayer(lakelist.filterMetadata('Hylak_id','equals', tmpid).style(styling), null, 'Lake polygon')

Answer 1

There's a bunch of bad practices in this script, but the main problem you're running into is that your first reduceRegion is using bestEffort=true, with a scale=30, and a Bbox that spans the entire planet. So you're asking for a scale that fits 1,784,452,180,388 into the default 10,000,000 pixel limit. That's a scale of 30m/pixel * (1784452180388/10000000) = 5,353,356m/pixel. At that scale, your polygon contains 0 pixels.

General Fixes:

• Don't ever use bestEffort; you have no idea what scale it's using.
• Don't use a global bounding box.
• Don't use tileScale unless you actually need it; you're asking for 256x more tiles which is just slowing things down.