在 Earth Engine/GIS 中提取样带数据

Question

我正在尝试研究如何在 Google 地球引擎中沿样带提取数据 - 最好是在定义的时间间隔内。在Jupyter notebook 中使用Python API 我已经映射了一些数据，缓冲了一个点（以定义感兴趣的区域），并映射了一个样带。我不知道我是否应该使用 ee 方法来提取沿lineString 约束的数据（我认为这不是一个形状？），或者我是否朝着错误的方向前进并且应该导出缓冲区域作为 GeoTIFF 来处理QGIS/ArcGIS 中的样带。

import ee
ee.Initialize()

def maskS2clouds(image):
  qa = image.select('QA60')

  # Bits 10 and 11 are clouds and cirrus, respectively.
  cloudBitMask = 1 << 10
  cirrusBitMask = 1 << 11

  # Both flags should be set to zero, indicating clear conditions.
  mask = qa.bitwiseAnd(cloudBitMask).eq(0) \
      .And(qa.bitwiseAnd(cirrusBitMask).eq(0))

  return image.updateMask(mask).divide(10000)

centre = ee.Geometry.Point(-115.435272, 35.584001,)
centreBuffer = centre.buffer(5000)

dataset = ee.ImageCollection('COPERNICUS/S2_SR') \
                  .filterDate('2020-07-01', '2020-07-31') \
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20)) \
                  .map(maskS2clouds)

visualization = {
  'min': 0.0,
  'max': 0.3,
  'bands': ['B4', 'B3', 'B2'], # Red, Green, BLue.  10m resolution.
}

Map.setCenter(-115.435272, 35.584001, 12)
# Map = geemap.Map(center=[35.584001, -115.435272], zoom=14)
Map.addLayer(dataset.mean(), visualization, 'RGB')
Map.addLayer(centre,
             {'color': 'black'},
             'Geometry [black]: point');
Map.addLayer(centreBuffer,
             {'color': 'red'},
             'Result [red]: point.buffer')

transect = ee.Geometry.LineString([[-115.4, 35.584001], [-115.45, 35.584001]]);
Map.addLayer(transect, {'color': 'Green'}, 'transect');

Map

Answer 1

您可能会发现吴秋生的这些资源很有帮助。

• https://www.youtube.com/watch?v=0TNXSs6fwrg
• https://geemap.org/common/#geemap.common.extract_transect
• https://geemap.org/notebooks/73_transect/

他的Common module API 有函数extract_transect(image, line, reducer='mean', n_segments=100, dist_interval=None, scale=None, crs=None, crsTransform=None, tileScale=1.0, to_pandas=False, **kwargs)

这会从图像中提取横断面。

从上面的第三个链接，Jupyter notebook 是：

import ee
import geemap
from bqplot import pyplot as plt

Map = geemap.Map()
Map.add_basemap("TERRAIN")
Map

image = ee.Image('USGS/SRTMGL1_003')
vis_params = {
  'min': 0,
  'max': 4000,
  'palette': ['006633', 'E5FFCC', '662A00', 'D8D8D8', 'F5F5F5']}
Map.addLayer(image, vis_params, 'SRTM DEM', True, 0.5)

# Use the drawing tool to draw any line on the map.
line = Map.user_roi
if line is None:
    line = ee.Geometry.LineString([[-120.223279, 36.314849], [-118.926969, 36.712192], [-117.202217, 36.756215]])
    Map.addLayer(line, {}, "ROI")   
Map.centerObject(line)

line.getInfo()

reducer='mean'  # Any ee.Reducer, e.g., mean, median, min, max, stdDev
transect = geemap.extract_transect(image, line, n_segments=100, reducer=reducer, to_pandas=True)

transect

fig = plt.figure()
plt.plot(transect['distance'], transect[reducer])
plt.xlabel('Distance')
plt.ylabel("Elevation")
plt.show()

JavaScript 版本是：

/***
 * Reduces image values along the given line string geometry using given reducer.
 * 
 * Samples image values using image native scale, or opt_scale
 */
function reduceImageProfile(image, line, reducer, scale, crs) {
  var length = line.length();
  var distances = ee.List.sequence(0, length, scale)
  var lines = line.cutLines(distances, ee.Number(scale).divide(5)).geometries();
  lines = lines.zip(distances).map(function(l) { 
    l = ee.List(l)
    
    var geom = ee.Geometry(l.get(0))
    var distance = ee.Number(l.get(1))
    
    geom = ee.Geometry.LineString(geom.coordinates())
    
    return ee.Feature(geom, {distance: distance})
  })
  lines = ee.FeatureCollection(lines)

  // reduce image for every segment
  var values = image.reduceRegions( {
    collection: ee.FeatureCollection(lines), 
    reducer: reducer, 
    scale: scale, 
    crs: crs
  })
  
  return values
}


// Define a line across the Olympic Peninsula, USA.

// Import a digital surface model and add latitude and longitude bands.
var elevImg = ee.Image('JAXA/ALOS/AW3D30/V2_2').select('AVE_DSM');

var profile = reduceImageProfile(elevImg, transect, ee.Reducer.mean(), 100)

print(ui.Chart.feature.byFeature(profile, 'distance', ['mean']))