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>嗯,用户想了解如何用Google Earth Engine(GEE)和随机森林算法进行水质反演。首先,我需要回顾一下GEE的基本功能和随机森林在遥感中的应用。GEE提供了大量的遥感数据集和处理能力,适合处理大范围的水质监测。随机森林作为一种机器学习算法,适合处理多源数据和非线性关系,常用于反演水质参数,比如叶绿素浓度、悬浮物含量等。 接下来,我应该整理实现步骤。通常,水质反演流程包括数据准备、特征选择、模型训练反演与验证。用户可能需要具体的代码示例,所以得准备一个结构清晰的代码框架,包括导入数据集、预处理、样本点提取模型训练和应用等步骤。 然后,注意事项部分需要考虑GEE的数据处理特点,比如影像的掩膜处理、波段选择,以及随机森林的参数调整。引用中的信息也很重要,比如用户提到的updateMask导致的背景问题,可能需要用栅格计算器处理,这点需要提及并引用相关文献。 用户可能对如何获取实测数据、特征工程优化、模型验证指标以及与其他算法的对比感兴趣,所以相关问题需要涵盖这些方面。同时,要确保回答中的LaTeX格式正确,比如变量用$...$,独立公式用$$...$$,并正确添加引用标识。 最后,要检查回答是否符合系统级指令,特别是结构清晰、步骤明确,并且引用恰当。例如,在代码示例后引用用户提供的案例,以及在处理背景问题时引用第二个引用源。确保所有技术细节准确无误,并且语言简洁易懂,帮助用户逐步理解和实现水质反演的过程。</think>### 如何使用Google Earth Engine (GEE) 和随机森林算法进行水质反演? #### **实现步骤** 1. **数据准备与预处理** - **遥感数据**:选择与水质参数(如叶绿素a、悬浮物、透明度等)相关的卫星影像(如Landsat-8、Sentinel-2)。例如: ```javascript var imageCollection = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') .filterDate('2020-01-01', '2020-12-31') .filterBounds(geometry); // geometry为研究区域 ``` - **实测数据**:收集与影像时间匹配的地面实测水质数据,格式为包含经纬度和参数值的CSV文件,需上传至GEE资产[^1]。 - **特征工程**:提取影像的波段反射率(如蓝、绿、红、近红外)、衍生指数(如NDWI、NDCI)作为输入特征: ```javascript var addIndices = function(image) { var ndwi = image.normalizedDifference(['SR_B3', 'SR_B5']).rename('NDWI'); var ndci = image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDCI'); return image.addBands(ndwi).addBands(ndci); }; ``` 2. **样本点匹配与数据集构建** - 将实测点与影像像素对齐,使用`ee.Image.sampleRegions`提取特征值: ```javascript var samples = image.select(['SR_B2','SR_B3','NDWI','NDCI']) // 选择特征波段 .sampleRegions({ collection: points, // 实测点集合 properties: ['Chl_a'], // 目标变量(如叶绿素a) scale: 30 }); ``` 3. **随机森林模型训练** - 使用GEE的`ee.Classifier.smileRandomForest`接口: ```javascript var trainedModel = ee.Classifier.smileRandomForest(100) // 100棵树 .setOutputMode('REGRESSION') .train({ features: samples, classProperty: 'Chl_a', inputProperties: ['SR_B2','SR_B3','NDWI','NDCI'] }); ``` 4. **水质参数反演与验证** - 应用模型至整幅影像,生成水质参数空间分布图: ```javascript var result = image.classify(trainedModel).rename('Chl_a_predicted'); ``` - 验证方法:划分训练集/测试集,计算均方根误差(RMSE): $$ \text{RMSE} = \sqrt{\frac{1}{n}\sum_{i=1}^n (y_i - \hat{y}_i)^2} $$ #### **注意事项** 1. **数据一致性**:实测数据与影像时间需匹配(±3天内最佳),避免因时间差异导致误差。 2. **特征优化**:通过相关性分析或变量重要性排序(如随机森林的`explain`方法)筛选关键特征。 3. **掩膜处理**:使用`updateMask`去除云、阴影等干扰区域,导出后若存在背景值问题,可用ArcGIS栅格计算器执行`影像*1`消除冗余背景[^2]。 #### **代码示例** ```javascript // 数据准备 var image = ee.Image('LANDSAT/LC08/C02/T1_L2/LC08_123032_20200101') .clip(geometry); var points = ee.FeatureCollection('users/your_asset/water_quality_points'); // 特征计算 var imageWithIndices = addIndices(image); // 样本提取 var samples = imageWithIndices.select(['SR_B2','SR_B3','NDWI','NDCI']) .sampleRegions({collection: points, properties: ['Chl_a'], scale: 30}); // 模型训练与预测 var model = ee.Classifier.smileRandomForest(100).train(samples, 'Chl_a'); var predicted = imageWithIndices.classify(model); // 导出结果 Export.image.toDrive({ image: predicted, description: 'Chl_a_map', scale: 30, region: geometry }); ```
评论 11
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

打赏作者

吕波涛.

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

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

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

打赏作者

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

抵扣说明:

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

余额充值