Google Earth Engine （GEE）——reduceRegion函数降低分辨率中出现错误计算的reducer.min从0变成了1

此星光明

发布于 2024-02-02 14:05:13

1850

发布于 2024-02-02 14:05:13

问题：

我目前正试图用reduceRegion函数找到一个二进制频段的最小值，也就是说，我想知道这个频段是否有0值。这个波段的空间细节非常粗糙，所以我想在降低的分辨率下运行（50米而不是原来的10米），以使它更有效率。然而，在50米的分辨率下，往往找不到正确的最小值。在10米的分辨率下，它可以正常工作，但这太耗费计算了。

在此脚本中，我们尝试检测低 S1 反向散射的簇，并创建一个这些集群的缓冲 (fastDistanceTransform) 掩码。我们添加这个mask到原始 S1 图像。我们想添加一个属性来表示是否低是否观察到反向散射簇，在整个过程中使用 reduceRegion分辨率为 50m。这是行不通的；它只适用于较低的分辨率10m，但这计算量更大。

解决方案：

当您在 50m 处运行 reduceRegion 时，10m 处的原始像素将使用“均值”采样进行聚合。所以 25 个原始像素将被转换为 1 个像素及其平均值。因此，即使只有一个非零像素，您得到的 50m 像素也不会为 0。这里具体的含义就是我们分辨率变粗的时候，就会出现原来很小的像素本来是0，但是随着统计范围的扩大，周围像素值只要有一个为1，那么就不会出现统计值为0的情况。你可以在https://developers.google.com/earth-engine/scale上查看比例在 GEE 中的工作原理。

另外，这里正确的方法是使用导出。每当您看到某些东西超时时，将其导出将允许它运行更长时间并访问更大的资源池。你可以导出任何内容，包括数字。只需要将它包装在一个特征集合中。这里我们的做法就是将其运算量过大的东西不要通过展示的形式显示，可以通过下载的形式导出到Google硬盘中。我们把我们需要打印的东西导出的过程就是把器装入一个feature中，然后再把逐个的feature放入到矢量集合中。就行了。

函数：

ee.Kernel.square(radius, units, normalize, magnitude)

Generates a square-shaped boolean kernel.

Arguments:

radius (Float):

The radius of the kernel to generate.

units (String, default: "pixels"):

The system of measurement for the kernel ('pixels' or 'meters'). If the kernel is specified in meters, it will resize when the zoom-level is changed.

normalize (Boolean, default: true):

Normalize the kernel values to sum to 1.

magnitude (Float, default: 1):

Scale each value by this amount.

Returns: Kernel

ee.Feature(geometry, properties)

Features can be constructed from one of the following arguments plus an optional dictionary of properties:

An ee.Geometry.
A GeoJSON Geometry.
A GeoJSON Feature.
A computed object: reinterpreted as a geometry if properties are specified, and as a feature if they aren't.

Arguments:

在本次错误修复中我们使用的第一个参数是没有的，因为我们只需要导出我们所需要的表格就行，这里的第一个研究区设定为null，第二个参数设定我们要导出的属性，这个案例中是min最小值。

geometry (ComputedObject|Feature|Geometry|Object):

A geometry or feature.

properties (Object, optional):

A dictionary of metadata properties. If the first parameter is a Feature (instead of a geometry), this is unused.

Returns: Feature

错误代码：

var geometry = 
    /* color: #d63000 */
    /* shown: false */
    ee.Geometry.Point([15.502206898267596, 1.7847823530702267]),
    geometry2 = 
    /* color: #98ff00 */
    /* shown: false */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[15.990200854789757, 1.8033025508734506],
          [15.990200854789757, 1.5383703103304462],
          [16.305371142875696, 1.5383703103304462],
          [16.305371142875696, 1.8033025508734506]]], null, false);


// Load a Sentinel-1 ImageCollection, take just one image to anayze
var s1_image = ee.ImageCollection('COPERNICUS/S1_GRD')
                    .filterDate('2020-03-01', '2020-05-01')
                    .filterBounds(geometry).first().clip(geometry2).select('VH')
                    
Map.addLayer(s1_image ,{min:-25,max:0}, 'S1 image')

// Define the pixels where a low backscatter is observed
var low_backscatter = s1_image.lt(-25)

// Get the CRS
var crs = s1_image.select(0).projection().crs()

// Identify the areas with large clusters of low backscatter pixels
var low_backscatter_clusters = low_backscatter.reduceNeighborhood({reducer:ee.Reducer.mean(), kernel:ee.Kernel.square(25), optimization:'boxcar'})
low_backscatter_clusters = low_backscatter_clusters.gt(0.5).rename('low_backscatter_clusters')

Map.addLayer(low_backscatter_clusters.reproject(crs, null, 10),{min:0,max:1}, 'Low backscatter clusters')


// Buffer the low backscatter clusters
var fastdistancetransform = low_backscatter_clusters.fastDistanceTransform(1000) //Note: distance is squared
var low_backscatter_clusters_buffered = fastdistancetransform.gt(Math.pow(1000/10, 2)).rename('low_backscatter_clusters') 

// Mask to the original image extent
low_backscatter_clusters_buffered = low_backscatter_clusters_buffered.updateMask(s1_image.select(0).abs())


Map.addLayer(low_backscatter_clusters_buffered.reproject(crs, null, 10), {min:0,max:1}, 'Low backscatter clusters buffered')


// Check if low backscatter clusters were detected in the image
var min = low_backscatter_clusters_buffered.reduceRegion({
  reducer: ee.Reducer.min(),
  geometry: s1_image.geometry(),
  scale: 50, //10  //It works at 10m resolution, but not at 50m resolution
  maxPixels: 1e13
}).get('low_backscatter_clusters')

var check = ee.Algorithms.If(ee.Number(min).eq(0), 'Yes', 'No')

// Add this information to the properties of the image
print('Did we observe low backscatter clusters?', check)
s1_image = s1_image.addBands(low_backscatter_clusters_buffered).set('Low_backscatter_clusters', check)

10和20分辨率都可以

但是超过30米分辨率的情况下就不行了

正确代码：

var geometry = 
    /* color: #d63000 */
    /* shown: false */
    ee.Geometry.Point([15.502206898267596, 1.7847823530702267]),
    geometry2 = 
    /* color: #98ff00 */
    /* shown: false */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[15.990200854789757, 1.8033025508734506],
          [15.990200854789757, 1.5383703103304462],
          [16.305371142875696, 1.5383703103304462],
          [16.305371142875696, 1.8033025508734506]]], null, false);

// Load a Sentinel-1 ImageCollection, take just one image to anayze
var s1_image = ee.ImageCollection('COPERNICUS/S1_GRD')
                    .filterDate('2020-03-01', '2020-05-01')
                    .filterBounds(geometry).first().clip(geometry2).select('VH')
                    
Map.addLayer(s1_image ,{min:-25,max:0}, 'S1 image')

// Define the pixels where a low backscatter is observed
var low_backscatter = s1_image.lt(-25)

// Get the CRS
var crs = s1_image.select(0).projection().crs()

// Identify the areas with large clusters of low backscatter pixels
var low_backscatter_clusters = low_backscatter.reduceNeighborhood({reducer:ee.Reducer.mean(), kernel:ee.Kernel.square(25), optimization:'boxcar'})
low_backscatter_clusters = low_backscatter_clusters.gt(0.5).rename('low_backscatter_clusters')

Map.addLayer(low_backscatter_clusters.reproject(crs, null, 10),{min:0,max:1}, 'Low backscatter clusters')


// Buffer the low backscatter clusters
var fastdistancetransform = low_backscatter_clusters.fastDistanceTransform(1000) //Note: distance is squared
var low_backscatter_clusters_buffered = fastdistancetransform.gt(Math.pow(1000/10, 2)).rename('low_backscatter_clusters') 

// Mask to the original image extent
low_backscatter_clusters_buffered = low_backscatter_clusters_buffered.updateMask(s1_image.select(0).abs())


Map.addLayer(low_backscatter_clusters_buffered.reproject(crs, null, 10), {min:0,max:1}, 'Low backscatter clusters buffered')


// Check if low backscatter clusters were detected in the image
var min = low_backscatter_clusters_buffered.reduceRegion({
  reducer: ee.Reducer.min(),
  geometry: s1_image.geometry(),
  scale: 10, //10  //It works at 10m resolution, but not at 50m resolution
  maxPixels: 1e13
}).get('low_backscatter_clusters')

var fc = ee.FeatureCollection([ee.Feature(null, {'min': min})]);

Export.table.toDrive({
  collection: fc, 
  description: 'min_value_export',
  folder: 'earthengine',
  fileNamePrefix: 'min',
  fileFormat: 'CSV'})