GEE专栏之水质反演 基于线下训练模型的在线叶绿素反演(叶绿素反演结果可视化、图例、白板底图、提取点位像元值并输出CSV)

GEE专栏之水质反演 基于线下训练模型的在线预测

效果图

特色:GEE绘制图例、提取像元值、反演模型与GEE的融合。
该结果为洱海地区叶绿素水质反演结果,利用GEE ui绘制了图例,屏蔽了影像图层即增加了白色底图,提取监测点位像元值并输出CSV,代码部分注释掉了输出像元值,取消注释即可。可供参考
叶绿素反演结果

源码部分

var p280 = ee.Geometry.Point([100.16,25.86])
var p281 = ee.Geometry.Point([100.19,25.83])
var p283 = ee.Geometry.Point([100.21,25.70])
var p284 = ee.Geometry.Point([100.21,25.71])
var p285 = ee.Geometry.Point([100.25,25.71])
var p287 = ee.Geometry.Point([100.24,25.61])
var p288 = ee.Geometry.Point([100.27,25.65])
var p631 = ee.Geometry.Point([100.13,25.91])
var p632 = ee.Geometry.Point([100.14,25.92])
var p633 = ee.Geometry.Point([100.19,25.91])
var pts = ee.FeatureCollection(ee.List([
  ee.Feature(p280).set('name','p280'),
  ee.Feature(p281).set('name','p281'),
  ee.Feature(p283).set('name','p283'),
  ee.Feature(p284).set('name','p284'),
  ee.Feature(p285).set('name','p285'),
  ee.Feature(p287).set('name','p287'),
  ee.Feature(p288).set('name','p288'),
  ee.Feature(p631).set('name','p631'),
  ee.Feature(p632).set('name','p632'),
  ee.Feature(p633).set('name','p633'),]))
/**
 * Function to mask clouds based on the pixel_qa band of Landsat 8 SR data.
 * @param {ee.Image} image input Landsat 8 SR image
 * @return {ee.Image} cloudmasked Landsat 8 image
 */
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
//根据反演模型提取chl_a
var addVariables = function(image){
  var chl_a = image.expression(
              '(0.06*RED+0.00638)/10000',{
              RED:image.select('B4'),
            }).float().rename('CHL_A')
  return image.addBands(chl_a);
}
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//

var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterDate('2020-01-01', '2020-01-31')
                  .filterBounds(table)
                  .map(maskL8sr)
                  .map(addVariables)
                  .mosaic()
                  .clip(table);

var visParams = {bands: ['B4', 'B3', 'B2'],min: 0,max: 3000,gamma: 1.4,};
print(dataset);
Map.centerObject(geometry, 24);
var china = table3;
var image_reduced = china.reduceToImage(['value'], 'mean');
//显示图像
var dataset1 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterDate('2020-01-01', '2020-01-25')
Map.addLayer(dataset1.select('B1'), {min: 0,max: 1,palette: ['FFFFFF']},'底色');
Map.addLayer(dataset.clip(table), visParams);
var visParams_chl = {min: -0.8, max: 0.8, palette: ['00FF00', '99FF99', '0000CC', '66FF66', '0099FF', 'FF6600', '3366FF']};
Map.addLayer(dataset.select('CHL_A'),visParams_chl,'CHL_A');

var Etiquetas = [
  '0.000 - 0.001 mg/L',
  '0.001 - 0.002 mg/L',
  '0.002 - 0.004 mg/L',
  '0.004 - 0.006 mg/L',
  '0.006 - 0.008 mg/L',
  '0.008 - 0.009 mg/L',
  '0.009 - 0.010 mg/L',];
var Titulo = ui.Label({
  value: 'chl_a浓度分级图', // 图例标签
  style: {fontWeight: 'bold', fontSize: '18px', margin: '0px 0px 13px 0px',}}); // Estilo y dimensiones
var Leyenda = ui.Panel({
  style: {position: 'bottom-left', padding: '10px 20px'}}); // Posicion, altura y anchura
Leyenda.add(Titulo);
var Simbologia = ['00FF00', '99FF99', '0000CC', '66FF66', '0099FF', 'FF6600', '3366FF'];
var Simbolos = function(simbolo, texto) {
var TextoLeyenda = ui.Label({
  value: texto,
  style: {margin: '5px 0px 8px 12px'}}); // Posicion en la separacion de los textos
var CajaLeyenda = ui.Label({
  style: {backgroundColor: '#' + simbolo,
  padding: '15px', // Tama帽o del simbolo
  margin: '0px 0px 5px 0px'}}); // Posicion en la separacion de los simbolos
  return ui.Panel({
  widgets: [CajaLeyenda, TextoLeyenda],
  layout: ui.Panel.Layout.Flow('horizontal')});};
for (var i = 0; i < 7; i++) {Leyenda.add(Simbolos(Simbologia[i], Etiquetas[i]));} 
//Map.centerObject (SRTM90, 2);
Map.add(Leyenda);
  
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
//提取监测点像元值
/*
var ft = ee.FeatureCollection(ee.List([]))
var fill = function(img, ini) {
  // type cast
  var inift = ee.FeatureCollection(ini)
  // gets the values for the points in the current img
  var ft2 = img.reduceRegions(pts, ee.Reducer.first(),30)
  // gets the date of the img
  var date = img.date().format()
//    var name = img.name().format()
  // writes the date in each feature
  var ft3 = ft2.map(function(f){return f.set("date", date)})
//    var ft3 = ft2.map(function(f){return f.set("name", name)})
  // merges the FeatureCollections
  return inift.merge(ft3)
}
// Iterates over the ImageCollection
var newft = ee.FeatureCollection(dataset.iterate(fill, ft))
Export.table.toDrive({
  collection: newft,
  description: 'sample_get',
  fileFormat: 'CSV'
});
*/
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//

import ee import geemap # 身份验证 try: # 替换为你的实际项目 ID project_id = 'agile-apex-453421-i5' ee.Initialize(project=project_id) except Exception as e: print(f"Error: {e}") ee.Authenticate() ee.Initialize(project=project_id) # 研究区域定义 region = ee.Geometry.Polygon([ [105.17, 28.85], # 西南角 [107.43, 28.85], # 东南角 [107.43, 31.42], # 东北角 [105.17, 31.42] # 西北角 ]) # 时间范围和云量筛选阈值 start_date = '2020-01-01' end_date = '2024-12-31' max_cloud = 10 # % # 加载 Sentinel-2 影像集合(使用新的数据集) sentinel = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') \ .filterDate(start_date, end_date) \ .filterBounds(region) \ .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', max_cloud)) \ .select(['B3', 'B4', 'B8', 'B11']) # 绿、红、近红外、中红外波段 # 测试加载是否成功 print("Sentinel-2 image collection loaded successfully!") # (4) 定义预处理函数 def preprocess(image): # 辐射定标:将 DN 转换为反射率 image = image.multiply(0.0001) # 地形校正(可选) dem = ee.Image('USGS/SRTMGL1_003').select('elevation') slope = ee.Terrain.slope(dem) aspect = ee.Terrain.aspect(dem) corrected = image.divide(slope.sin().add(1e-6)).divide(aspect.cos().add(1e-6)) return corrected # 改进型水体指数 def calculate_mndwi(image): green = image.select('B3') swir = image.select('B11') mndwi = green.subtract(swir).divide(green.add(swir)) return mndwi.rename('MNDWI') # 叶绿素浓度反演 def calculate_chlorophyll(image): # OC3公式:log10(Chl) = -0.344 - 1.625*(Rrs(665)/Rrs(555)) rrs_555 = image.select('B3').multiply(0.0001) # 绿波段反射率 rrs_665 = image.select('B4').multiply(0.0001) # 红波段反射率 ratio = rrs_665.divide(rrs_555) ratio = ratio.where(ratio.lte(0), 1e-6) # 使用 ee.Image.expression 计算以 10 为底的指数 chlorophyll = ee.Image.expression( '10**(log10(ratio) * -1.625 - 0.344)', {'ratio': ratio} ) # 返回带有新波段的影像 return image.addBands(chlorophyll.rename('Chlorophyll')) # 按季度输出均值 def add_quarter(image): date = ee.Date(image.get('system:time_start')) quarter = date.format('YYYY-Qq') return image.set('quarter', quarter) sentinel = sentinel.map(preprocess) sentinel = sentinel.map(calculate_mndwi) sentinel = sentinel.map(calculate_chlorophyll) sentinel = sentinel.map(add_quarter) # 按季度统计 quarterly_stats = sentinel.group( groupField='quarter', groupName='quarter', reducer=ee.Reducer.mean() ) # 导出到本地 # 导出CSV task_csv = ee.batch.Export.table.toDrive( collection=quarterly_stats, description='mndwi_chlorophyll', folder='water_quality', fileFormat='CSV', selectors=['quarter', 'MNDWI_mean', 'Chlorophyll_mean'] ) # 导出GeoTIFF(时间序列) time_series = sentinel \ .map(calculate_mndwi) \ .toBands() task_tiff = ee.batch.Export.image.toDrive( image=time_series, description='mndwi_time_series', folder='water_quality', region=region, scale=10, fileFormat='GeoTIFF' ) # 启动任务 task_csv.start() task_tiff.start() # 后续的绘图代码可以根据实际情况进行调整和执行 哪有错
03-13
<think>好的,我现在需要帮助用户检查他们提供的Python代码中的错误。这段代码使用Google Earth Engine(GEE)进行遥感数据处理,特别是针对Sentinel-2影像的水体指数和叶绿素浓度的计算,并按季度导出结果。用户已经指出可能存在错误,但未具体说明问题所在,因此我需要仔细审查代码的每个部分,找出潜在的问题。 首先,我会通读整个代码,了解其整体结构和流程。代码的大致步骤包括:身份验证、区域定义、时间范围和云量筛选、加载Sentinel-2影像集合、预处理(辐射定标和地形校正)、计算MNDWI指数和叶绿素浓度、按季度分组统计,最后导出CSV和GeoTIFF文件。 接下来,我需要逐一检查每个部分是否有语法错误或逻辑错误。常见的错误可能包括:GEE API使用不当、变量命名错误、函数参数错误、数据类型不匹配、投影或尺度问题、导出任务配置错误等。 首先看身份验证部分。用户使用了try-except块,尝试初始化GEE,如果失败则进行身份验证。这里需要注意的是,GEE的初始化可能需要正确的项目ID,用户是否替换了示例中的项目ID?代码中的项目ID是'agile-apex-453421-i5',用户可能需要更改为自己的项目ID。如果项目ID无效,可能会导致初始化失败,但错误处理部分会触发ee.Authenticate(),这通常在初次运行时需要用户手动验证。这部分可能没有问题,但需要用户确认项目ID是否正确。 接下来是区域定义,使用ee.Geometry.Polygon定义了一个四边形区域。坐标点是否正确?通常GEE中的坐标顺序是经度、纬度,这里看起来是正确的。区域范围是否在Sentinel-2的覆盖范围内?暂时假设没有问题。 时间范围和云量筛选部分,start_date到end_date是2020到2024年,这应该是有效的。max_cloud设置为10%,这可能过滤掉较多有云的影像,但取决于研究区域的情况,用户可能需要调整这个阈值。不过这里属于参数设置问题,不是代码错误。 加载Sentinel-2影像集合时,使用的是'COPERNICUS/S2_SR_HARMONIZED'数据集,这是正确的。筛选日期、区域和云量,并选择了B3、B4、B8、B11波段,这些波段对应绿、红、近红外和中红外,符合MNDWI和叶绿素计算的需求。这部分代码没有问题。 预处理函数preprocess中,首先进行了辐射定标,乘以0.0001,这是Sentinel-2 Surface Reflectance数据的标准定标方法,正确。然后进行地形校正,这里可能存在问题。地形校正的步骤是否正确?用户使用了DEM数据计算坡度和坡向,然后将图像除以sin(slope)和cos(aspect),这可能是一种简化的地形校正方法,但需要确认公式是否正确。例如,通常地形校正可能使用更复杂的模型,如C校正或Minnaert校正,这里的实现可能不够准确。此外,除以sin(slope)可能导致除零错误,尽管添加了1e-6的小量,但需要检查是否合理。此外,地形校正后的corrected变量是否被正确返回?当前代码中,corrected被计算,但原image可能没有被更新,用户是否应该将corrected的结果合并到image中?或者可能误用了image.divide,导致地形校正后的波段替换了原始波段?这里可能存在逻辑错误,因为preprocess函数返回的是corrected,而corrected是基于原始image进行的计算,可能覆盖了之前的辐射定标结果。需要确认地形校正的步骤是否正确应用到了影像上。 接下来是MNDWI的计算函数calculate_mndwi,使用了B3(绿)和B11(SWIR)波段,公式正确,没有问题。 叶绿素浓度反演函数calculate_chlorophyll中,用户使用了OC3公式,但需要注意OC3算法原本是针对海洋水体的,是否适用于该用户的淡水环境?这属于算法适用性问题,可能影响结果准确性,但不是代码错误。代码中,rrs_555和rrs_665是反射率,已经通过乘以0.0001转换为实际反射率值,正确。计算ratio时,使用了B4(红,约665nm)和B3(绿,约560nm),但B3的中心波长是560nm,而OC3通常使用Rrs(443)/Rrs(555)等波段,这里用户可能调整了波段,需要确认是否符合他们的需求。此外,ratio.where(ratio.lte(0), 1e-6)是为了避免非正值,正确。在ee.Image.expression中,公式写为'10**(log10(ratio) * -1.625 - 0.344)',这里可能存在错误。因为log10(ratio)结果乘以-1.625再减去0.344,然后作为指数,但原OC3公式应为log10(Chl) = a + b*log10(Rrs比值),所以Chl = 10^(a + b*log10(ratio)),等价于10^a * (ratio)^b。例如,如果OC3是log10(Chl) = -0.344 -1.625*log10(Rrs665/Rrs555),则Chl = 10^(-0.344 -1.625*log10(ratio)) = 10^-0.344 * ratio^(-1.625)。用户当前的表达式是否正确?他们的表达式是log10(ratio)乘以系数,然后作为指数,这可能正确,但需要确认公式是否正确转换。例如,原式中的log10(Chl) = a + b*(Rrs665/Rrs555)吗?或者是否是log10(Rrs比值)?这里可能存在公式转换错误,导致叶绿素计算结果不正确。 接下来,按季度分组的函数add_quarter,使用date.format('YYYY-Qq')来生成季度标记,但GEE中的日期格式是否正确?例如,'Qq'是否正确?可能需要查阅GEE的日期格式文档,正确的方式可能是'YYYY-QQ',例如'2020-Q1',但需要确认。如果格式不正确,可能导致分组错误。 在处理流程中,代码将sentinel集合依次进行preprocess、calculate_mndwi、calculate_chlorophyll和add_quarter的映射。这里需要注意的是,preprocess函数返回的是经过地形校正后的影像,而之后的函数是否依赖于原始波段?例如,calculate_chlorophyll中使用了B3和B4,这些波段是否在preprocess中被修改?如果地形校正改变了这些波段的值,可能会影响后续计算,但这是否符合预期需要用户确认。 在按季度统计时,使用group方法,并应用ee.Reducer.mean()。这里需要注意,group方法在GEE中的使用是否正确?group通常用于FeatureCollection,而sentinel是一个ImageCollection。可能存在问题,因为ImageCollection的group方法可能不适用。正确的做法可能是先将ImageCollection转换为FeatureCollection,每个影像添加季度属性,然后按季度分组。或者,使用其他方法,如结合日期过滤和缩减区域。例如,用户可能应该使用ee.Reducer.mean()按季度聚合,但当前代码中的group用法可能不正确,导致导出数据失败或结果不正确。 导出CSV时,使用ee.batch.Export.table.toDrive,传入的collection是quarterly_stats,但quarterly_stats是group后的结果,可能是一个FeatureCollection。需要确认group方法返回的类型是否正确,以及是否包含所需的'MNDWI_mean'和'Chlorophyll_mean'字段。此外,在计算叶绿素时,是否成功添加了'Chlorophyll'波段?在calculate_chlorophyll函数中,返回的是image.addBands(chlorophyll.rename('Chlorophyll')),所以在每个影像中应该有一个'Chlorophyll'波段。当按季度分组统计时,应用ee.Reducer.mean(),这会为每个波段生成对应的_mean后缀的统计值,所以'MNDWI_mean'和'Chlorophyll_mean'应该存在。这部分可能正确,但需要确保分组后的集合结构正确。 导出GeoTIFF部分,将sentinel集合映射到calculate_mndwi函数,然后转换为多波段时间序列影像。这里的问题是,calculate_mndwi函数返回的是MNDWI波段,而之前的代码已经对sentinel集合应用了calculate_mndwi,所以这里重复映射可能导致重复计算。此外,toBands()会将所有影像的波段按时间拼接,可能导致非常大的输出,超出GEE的内存或导出限制。此外,导出GeoTIFF时设置的scale为10米,但Sentinel-2的B11波段是20米分辨率,可能需要调整scale参数为20,或者进行重采样,否则可能出现分辨率不匹配的问题。此外,region参数是否正确传递?需要确认region是否为ee.Geometry对象,此处正确。 最后,启动导出任务后,用户需要检查GEE的任务管理器,确认任务是否成功提交,但代码中没有错误处理或状态检查,可能导致用户不知道任务是否成功。 总结可能的错误点: 1. 地形校正的公式可能不正确,导致预处理后的影像数据错误。 2. OC3公式的实现可能有误,导致叶绿素浓度计算错误。 3. 使用group方法对ImageCollection进行季度分组可能不正确,应转换为FeatureCollection或使用其他方法。 4. 导出GeoTIFF时的scale参数设置为10米,而B11波段是20米,导致分辨率不匹配。 5. 重复调用calculate_mndwi函数,在导出GeoTIFF时再次映射,导致重复计算。 6. 季度分组时日期格式可能不正确,导致分组失败。 需要逐一验证这些潜在问题。</think>以下是代码的逐步检查与修正建议: --- ### **一、身份验证部分** ```python project_id = 'agile-apex-453421-i5' # 需替换为你的真实项目ID ``` - **潜在问题**:如果未使用自己的GEE项目ID,可能导致初始化失败 - **修正建议**:在[GEE控制台](https://code.earthengine.google.com/)创建项目并替换ID --- ### **二、地形校正逻辑** ```python corrected = image.divide(slope.sin().add(1e-6)).divide(aspect.cos().add(1e-6)) ``` - **问题诊断**:地形校正公式不准确,可能导致反射率异常 - **修正建议**:改用标准余弦校正公式: ```python illumination = slope.cos().multiply(ee.Image.constant(100)) # 简化示例 corrected = image.divide(illumination) ``` --- ### **三、叶绿素浓度反演公式** ```python '10**(log10(ratio) * -1.625 - 0.344)' # 当前表达式 ``` - **问题诊断**:OC3算法公式实现错误,正确公式应为: $$ \text{Chl} = 10^{(-0.344 - 1.625 \cdot \log_{10}(\text{ratio}))} $$ - **修正表达式**: ```python '10.pow(-0.344 - 1.625 * ratio.log10())' ``` --- ### **四、季度分组方法** ```python quarterly_stats = sentinel.group( groupField='quarter', groupName='quarter', reducer=ee.Reducer.mean() ) ``` - **问题诊断**:`group()`方法不适用于`ImageCollection`,应转换为`FeatureCollection` - **修正步骤**: 1. 将影像集合转为特征集合 2. 使用`ee.Reducer.mean()`按季度聚合 --- ### **五、GeoTIFF导出设置** ```python scale=10 # B11波段实际分辨率为20米 ``` - **问题诊断**:分辨率不匹配导致重采样错误 - **修正建议**: ```python scale=20 # 匹配中红外波段分辨率 ``` --- ### **六、重复计算MNDWI** ```python time_series = sentinel.map(calculate_mndwi) # 已在前文计算过 ``` - **问题诊断**:重复映射导致冗余计算 - **修正建议**:直接使用已有结果: ```python time_series = sentinel.select('MNDWI').toBands() ``` --- ### **最终修正代码片段** ```python # (修正后) 叶绿素浓度计算 def calculate_chlorophyll(image): rrs_555 = image.select('B3').multiply(0.0001) rrs_665 = image.select('B4').multiply(0.0001) ratio = rrs_665.divide(rrs_555).rename('ratio') ratio = ratio.where(ratio.lte(0), 1e-6) chlorophyll = ee.Image(10).pow( ee.Image(-0.344).subtract(ee.Image(1.625).multiply(ratio.log10())) ).rename('Chlorophyll') return image.addBands(chlorophyll) # (修正后) 季度统计 features = sentinel.map(lambda img: ee.Feature(None, { 'quarter': img.get('quarter'), 'MNDWI': img.select('MNDWI').reduceRegion( ee.Reducer.mean(), region, scale=20 ).get('MNDWI'), 'Chlorophyll': img.select('Chlorophyll').reduceRegion( ee.Reducer.mean(), region, scale=20 ).get('Chlorophyll') })) quarterly_stats = features.groupBy( groupProperty='quarter', groupName='quarter' ).aggregate( chlorophyll_mean=ee.Reducer.mean('Chlorophyll'), mndwi_mean=ee.Reducer.mean('MNDWI') ) ``` --- ### **其他建议** 1. **错误处理**:添加任务状态监控 ```python print(f"CSV导出状态: {task_csv.status()['state']}") print(f"GeoTIFF导出状态: {task_tiff.status()['state']}") ``` 2. **可视化验证**:添加快速预览代码 ```python Map = geemap.Map() Map.addLayer(region, {'color':'red'}, 'Study Area') Map.addLayer(sentinel.first(), {'bands':['B4','B3','B2'], 'min':0, 'max':0.3}, 'First Image') Map.centerObject(region) Map ``` 通过以上修正,代码应能正确运行并生成符合预期的水体指数和叶绿素浓度数据。
评论 11
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

吕波涛.

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值