【GEE】時間序列多源遙感數據隨機森林回歸預測|反演|驗證|散點圖|完整代碼

實驗介紹

分類和回歸之間的主要區別在于,在分類中,我們的預測目標是離散的類別,而在回歸中,預測目標是連續的預測值。

本實驗的研究區域位于佛蒙特州的埃塞克斯郡,使用訓練數據來模擬土壤氧化還原深度,然后生成準確度評估圖表和統計數據。(數據僅供實驗使用,不代表真實值)

實驗目標

  1. 隨機森林回歸

  2. GEE 圖表繪制

實驗數據

VT_boundary.shp – shapefile 表示感興趣的示例區域

VT_pedons.shp – 佛蒙特州埃塞克斯縣訓練數據的 shapefile

var VT_boundary = ee.FeatureCollection("projects/ee-yelu/assets/VT_boundary"),VT_pedons = ee.FeatureCollection("projects/ee-yelu/assets/essex_pedons_all");

實驗環境

Chrome瀏覽器

earth engine賬號

目錄

第 1 部分:合成時間序列多參數影像數據

第 2 部分:準備訓練/驗證數據

第 3 部分:運行隨機森林回歸

第 4 部分:向地圖添加回歸,創建圖例

第 5 部分:創建模型評估統計數據和圖表

第 6 部分:驗證

第 7 部分:導出

第 8 部分:討論

時間序列Sentinel-1、Sentinel-2影像預處理

上傳矢量數據到earth engine
確保您已將VT_boundary.shp文件上傳到您的assets文件夾并將其導入到您的腳本中。確保將其從table重命名為VT_boundary,然后將代碼保存到本地存儲庫。

多時相Sentinel-2影像預處理
因為研究區域位于不同的地理區域,因此使用earth engine 加載自定義矢量時

需要準確地定義矢量文件的投影。

// 研究區
var roi = VT_boundary.union(); 
Map.addLayer(roi, {}, 'shp', false);
var crs = 'EPSG:4326'; // EPSG number for output projection. 32618 = WGS84/UTM Zone 18N. For more info- http://spatialreference.org/ref/epsg/

S2 影像去云與合成

function maskS2clouds(image) {var qa = image.select('QA60');// Bits 10 and 11 are clouds and cirrus, respectively.var cloudBitMask = 1 << 10;var cirrusBitMask = 1 << 11;// Both flags should be set to zero, indicating clear conditions.var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));return image.updateMask(mask).divide(10000);
}

計算多種植被指數作為特征

// NDVI
function NDVI(img) {var ndvi = img.expression("(NIR-R)/(NIR+R)", {"R": img.select("B4"),"NIR": img.select("B8"),});return img.addBands(ndvi.rename("NDVI"));
}// NDWI
function NDWI(img) {var ndwi = img.expression("(G-MIR)/(G+MIR)", {"G": img.select("B3"),"MIR": img.select("B8"),});return img.addBands(ndwi.rename("NDWI"));
}// NDBI
function NDBI(img) {var ndbi = img.expression("(SWIR-NIR)/(SWIR-NIR)", {"NIR": img.select("B8"),"SWIR": img.select("B12"),});return img.addBands(ndbi.rename("NDBI"));
}//SAVI
function SAVI(image) {var savi = image.expression('(NIR - RED) * (1 + 0.5)/(NIR + RED + 0.5)', {'NIR': image.select('B8'),'RED': image.select('B4')}).float();return image.addBands(savi.rename('SAVI'));
}//IBI 
function IBI(image) {// Add Index-Based Built-Up Index (IBI)var ibiA = image.expression('2 * SWIR1 / (SWIR1 + NIR)', {'SWIR1': image.select('B6'),'NIR': image.select('B5')}).rename(['IBI_A']);var ibiB = image.expression('(NIR / (NIR + RED)) + (GREEN / (GREEN + SWIR1))', {'NIR': image.select('B8'),'RED': image.select('B4'),'GREEN': image.select('B3'),'SWIR1': image.select('B11')}).rename(['IBI_B']);var ibiAB = ibiA.addBands(ibiB);var ibi = ibiAB.normalizedDifference(['IBI_A', 'IBI_B']);return image.addBands(ibi.rename(['IBI']));
}//RVI
function RVI(image) {var rvi = image.expression('NIR/Red', {'NIR': image.select('B8'),'Red': image.select('B4')});return image.addBands(rvi.rename('RVI'));
}//DVI
function DVI(image) {var dvi = image.expression('NIR - Red', {'NIR': image.select('B8'),'Red': image.select('B4')}).float();return image.addBands(dvi.rename('DVI'));
}

逐月合成Sentinel-2、Sentinel-1影像,并計算植被指數(通過for循環實現逐月合成,可以根據需要修改為自定義的時間序列)

// 創建var inStack = ee.Image()
for (var i = 0; i < 3; i += 1) {var start = ee.Date('2019-03-01').advance(30 * i, 'day');print(start)var end = start.advance(30, 'day');var dataset = ee.ImageCollection('COPERNICUS/S2_SR').filterDate(start, end).filterBounds(roi).filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',75)).map(maskS2clouds).map(NDVI).map(NDWI).map(NDBI).map(SAVI).map(IBI).map(RVI).map(DVI);var inStack_monthly = dataset.median().clip(roi);// 可視化var visualization = {min: 0.0,max: 0.3,bands: ['B4', 'B3', 'B2'],};Map.addLayer(inStack_monthly, visualization, 'S2_' + i, false);// // 紋理特征var B8 = inStack_monthly.select('B8').multiply(100).toInt16();var glcm = B8.glcmTexture({size: 3});var contrast = glcm.select('B8_contrast');var var_ = glcm.select('B8_var');var savg = glcm.select('B8_savg');var dvar = glcm.select('B8_dvar');inStack_monthly = inStack_monthly.addBands([contrast, var_, savg, dvar])// Sentinel-1 var imgVV = ee.ImageCollection('COPERNICUS/S1_GRD').filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV')).filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH')).filter(ee.Filter.eq('instrumentMode', 'IW')).filterBounds(roi).map(function(image) {var edge = image.lt(-30.0);var maskedImage = image.mask().and(edge.not());return image.updateMask(maskedImage);});var img_S1_asc = imgVV.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));var VV_img = img_S1_asc.filterDate(start, end).select("VV").median().clip(roi);var VH_img = img_S1_asc.filterDate(start, end).select("VH").median().clip(roi);inStack_monthly = inStack_monthly.addBands([VV_img, VH_img])print("img_S2_monthly", inStack_monthly)inStack = inStack.addBands(inStack_monthly.select(inStack_monthly.bandNames()))
}inStack = inStack.select(inStack.bandNames().remove("constant"))
print(inStack)

在控制臺上輸出,堆疊后的多波段image:

print("Predictor Layers:", inStack);
inStack.reproject(crs, null, 30);
Map.centerObject(inStack)

單擊運行并耐心等待(數據量比較大因此耗時較長)一個名為“ Predictor Layers ”的image對象將出現在控制臺中。單擊“ Image ”旁邊的向下箭頭,然后單擊“ bands ”旁邊的向下箭頭以檢查此對象。可以看出,我們創建了多時相多參數的遙感影像

準備訓練/驗證數據

A. 加載 AOI pedons shapefile

在開始之前,需要將樣本數據VT_pedons.shp作為assets加載到 GEE ,并導入到我們的代碼中,以便接下來在回歸中使用。

轉到左側面板中的assets選項卡,找到shapefile,然后單擊它,此時會彈出預覽:

圖片

導航到“features”選項卡并瀏覽 shapefile 的不同屬性。本實驗目標預測土壤氧化還原深度,因此“ REDOX_CM ”是我們需要預測的變量。

單擊右下角的藍色IMPORT按鈕將其添加到我們的代碼中。

此時會看到一個table變量已添加到頂部的導入列表中。我們將新 shapefile 的名稱從默認的table更改為VT_pedons。

現在您應該有兩個import:VT_boundary和VT_pedons。

B. 使用 pedons 準備訓練和驗證數據

“VT_pedons”是用于回歸的樣本數據,我們現在將其轉為ee.FeatureCollection()類型,以便我們在GEE中調用。

var trainingFeatureCollection = ee.FeatureCollection(VT_pedons, 'geometry');

接下來我們開始用隨機森林做回歸

運行隨機森林回歸

A. 選擇需要使用的波段

出于本練習的目的,我們剛剛堆疊了多個波段,然而在一些研究中,可能僅需要某些波段參與回歸,因此可以通過

ee.Image().select("波段名稱")

來篩選需要的波段。

這里,我們仍然使用全部的波段進行回歸分析。

var bands = inStack.bandNames(); //All bands on included here

B.提取訓練數據

我們的下一步是根據樣本點坐標,提取影像上相應位置的多個波段的數據。

var training = inStack.reduceRegions({collection:trainingFeatureCollection,reducer :ee.Reducer.mean(),scale:100,tileScale :5,crs: crs
});
print(training)
var trainging2 = training.filter(ee.Filter.notNull(bands));

注意:我們將訓練數據的空間分辨率(scale)設置為 30 m——使其與我們之前應用的預測層重投影保持一致。
如果有興趣進一步探索sampleRegions命令,只需在左側面板的Docs搜索欄中輸入“ ee.Image.sampleRegions ”即可。

然后加載訓練數據,將80%/20% 用于訓練/驗證

// 在training要素集中增加一個random屬性,值為0到1的隨機數
var withRandom = trainging2.randomColumn({columnName:'random',seed:2,distribution: 'uniform',
});var split = 0.8; 
var trainingData = withRandom.filter(ee.Filter.lt('random', split));
var validationData = withRandom.filter(ee.Filter.gte('random', split));

以上代碼為回歸添加了訓練和驗證數據,

C. 運行 RF 分類器

然后,我們使用訓練數據來創建隨機森林分類器。盡管我們執行的是回歸,而不是分類,這仍然被稱為classifier。

// train the RF classification model
var classifier = ee.Classifier.smileRandomForest(100, null, 1, 0.5, null, 0).setOutputMode('REGRESSION').train({features: trainingData, classProperty: 'REDOX_CM', inputProperties: bands});
print(classifier);

注意’ setOutputMode’一定設置為’ REGRESSION '的,該參數對于在 GEE 中運行不同類型的隨機森林模型至關重要。

對于隨機森林超參數的設置可以查看GEE Docs,描述如下:

圖片

最后,現在我們將使用剛剛創建的分類器對圖像進行分類。

// Classify the input imagery.
var regression = inStack.classify(classifier, 'predicted');
print(regression);

值得一提的事,在 Earth Engine 中,即使我們正在執行回歸任務,但仍然使用的事classify()方法。

到目前為止,我們已經創建了一個空間回歸模型,但我們還沒有將它添加到我們的地圖中,所以如果您運行此代碼,您的控制臺或地圖中不會出現任何新內容(記得順手ctrl+s)

向地圖添加回歸結果,創建圖例

A. 加載并定義一個連續調色板

由于我們的回歸是對連續變量進行分類,因此我們不需要像分類時那樣為每個類選擇顏色。相反,我們將使用預制的調色板——要訪問這些調色板,我們需要加載庫。為此,請導入此鏈接以將模塊添加到您的 Reader 存儲庫

var palettes = require('users/gena/packages:palettes');

如果您想在將來選擇不同的連續調色板,請訪問此 URL。

現在我們已經加載了所需的庫,我們可以為回歸輸出定義一個調色板。在“選擇并定義調色板”的注釋下,粘貼:

var palette = palettes.crameri.nuuk[25];

B. 將最終回歸結果添加到地圖

現在我們已經定義了調色板,我們可以將結果添加到地圖中。

// Display the input imagery and the regression classification.
// get dictionaries of min & max predicted value
var regressionMin = (regression.reduceRegion({reducer: ee.Reducer.min(),scale: 100,crs: crs,bestEffort: true,tileScale: 5,geometry: roi
}));
var regressionMax = (regression.reduceRegion({reducer: ee.Reducer.max(),scale: 100,crs: crs,bestEffort: true,tileScale: 5,geometry: roi
}));// Add to map
var viz = {palette: palette,min: regressionMin.getNumber('predicted').getInfo(),max: regressionMax.getNumber('predicted').getInfo()
};
Map.addLayer(regression, viz, 'Regression');

如您所見,在地圖上顯示回歸的結果比顯示分類要復雜一些。這是因為該代碼的第一部分是為我們的可視化計算適當的最小值和最大值——它只是查找和使用預測氧化還原深度的最高和最低值。在后續計算中,當您使用此代碼對不同的連續變量進行建模時,它會自動為您的可視化選擇合適的值。

tileScale參數調整用于計算用于在地圖上適當顯示回歸的最小值和最大值的內存。如果您在此練習中遇到內存問題時,可以增加此參數的值。您可以在 GEE 指南中了解有關使用tileScale 和調試的更多信息。

制作圖例,將其添加到地圖

圖片

在地圖上顯示圖例總是很有用的,尤其是在處理各種顏色時。

以下代碼可能看起來讓人頭大,但其中大部分只是創建圖例的結構和其他美化細節。

// Create the panel for the legend items.
var legend = ui.Panel({style: {position: 'bottom-left',padding: '8px 15px'}
});// Create and add the legend title.
var legendTitle = ui.Label({value: 'Legend',style: {fontWeight: 'bold',fontSize: '18px',margin: '0 0 4px 0',padding: '0'}
});
legend.add(legendTitle);// create the legend image
var lon = ee.Image.pixelLonLat().select('latitude');
var gradient = lon.multiply((viz.max - viz.min) / 100.0).add(viz.min);
var legendImage = gradient.visualize(viz);// create text on top of legend
var panel = ui.Panel({widgets: [ui.Label(viz['max'])],
});legend.add(panel);// create thumbnail from the image
var thumbnail = ui.Thumbnail({image: legendImage,params: {bbox: '0,0,10,100',dimensions: '10x200'},style: {padding: '1px',position: 'bottom-center'}
});// add the thumbnail to the legend
legend.add(thumbnail);// create text on top of legend
var panel = ui.Panel({widgets: [ui.Label(viz['min'])],
});legend.add(panel);
Map.add(legend);// Zoom to the regression on the map
Map.centerObject(roi, 11);

圖片

創建模型評估統計數據

可視化評估工具

數據可視化是評估模型性能的一個極其重要的方法,很多時候通過可視化的方式看結果會容易得多。

我們要制作的第一個圖是具有特征重要性的直方圖。這是一個有用的圖標,尤其是當我們在模型中使用多個個預測層時。它使我們能夠查看哪些變量對模型有幫助,哪些變量可能沒有。

// Get variable importance
var dict = classifier.explain();
print("Classifier information:", dict);
var variableImportance = ee.Feature(null, ee.Dictionary(dict).get('importance'));
// Make chart, print it
var chart =ui.Chart.feature.byProperty(variableImportance).setChartType('ColumnChart').setOptions({title: 'Random Forest Variable Importance',legend: {position: 'none'},hAxis: {title: 'Bands'},vAxis: {title: 'Importance'}});
print(chart);

運行代碼后,您應該有類似下圖所示的內容:

圖片

接下來,我們將制作一個直方圖,顯示在每個氧化還原深度預測我們研究區域中有多少個像素。這是評估預測值分布的有用視覺效果。

// Set chart options
var options = {lineWidth: 1,pointSize: 2,hAxis: {title: 'Redox (cm)'},vAxis: {title: 'Num of pixels'},title: 'Number of pixels by redox depth'
};
var regressionPixelChart = ui.Chart.image.histogram({image: ee.Image(regression),region: roi,scale:100}).setOptions(options);
print(regressionPixelChart);

圖片

最后,我們將制作的最后一個圖是預測值和真實值的散點圖。這些對于查看模型的擬合情況十分有幫助,因為它從回歸圖像(預測值)中獲取樣本點,并將其與訓練數據(真實值)進行對比。

// Get predicted regression points in same location as training data
var predictedTraining = regression.sampleRegions({collection: trainingData,scale:100,geometries: true,projection:crs
});
// Separate the observed (REDOX_CM) and predicted (regression) properties
var sampleTraining = predictedTraining.select(['REDOX_CM', 'predicted']);
// Create chart, print it
var chartTraining = ui.Chart.feature.byFeature({features:sampleTraining, xProperty:'REDOX_CM', yProperties:['predicted']}).setChartType('ScatterChart').setOptions({title: 'Predicted vs Observed - Training data ',hAxis: {'title': 'observed'},vAxis: {'title': 'predicted'},pointSize: 3,trendlines: {0: {showR2: true,visibleInLegend: true},1: {showR2: true,visibleInLegend: true}}});
print(chartTraining);

圖片

注意:如果您將鼠標懸停在該圖的右上角,您也可以看到 R^2 值。(R^2 值會顯示在圖上)

計算均方根誤差 (RMSE)

我們將計算RMSE來評估我們的回歸對訓練數據的影響。

// Compute RSME
// Get array of observation and prediction values 
var observationTraining = ee.Array(sampleTraining.aggregate_array('REDOX_CM'));
var predictionTraining = ee.Array(sampleTraining.aggregate_array('predicted'));
print('observationTraining', observationTraining)
print('predictionTraining', predictionTraining)// Compute residuals
var residualsTraining = observationTraining.subtract(predictionTraining);
print('residualsTraining', residualsTraining)
// Compute RMSE with equation, print it
var rmseTraining = residualsTraining.pow(2).reduce('mean', [0]).sqrt();
print('Training RMSE', rmseTraining);

但是,僅通過查看我們的訓練數據,我們無法真正了解我們的數據表現如何。我們將對我們的驗證數據執行類似的評估,以了解我們的模型在未用于訓練它的數據上的表現如何。

// Get predicted regression points in same location as validation data
var predictedValidation = regression.sampleRegions({collection: validationData,scale:100,geometries: true
});
// Separate the observed (REDOX_CM) and predicted (regression) properties
var sampleValidation = predictedValidation.select(['REDOX_CM', 'predicted']);
// Create chart, print it
var chartValidation = ui.Chart.feature.byFeature(sampleValidation, 'predicted', 'REDOX_CM').setChartType('ScatterChart').setOptions({title: 'Predicted vs Observed - Validation data',hAxis: {'title': 'predicted'},vAxis: {'title': 'observed'},pointSize: 3,trendlines: {0: {showR2: true,visibleInLegend: true},1: {showR2: true,visibleInLegend: true}}});
print(chartValidation);// Get array of observation and prediction values 
var observationValidation = ee.Array(sampleValidation.aggregate_array('REDOX_CM'));
var predictionValidation = ee.Array(sampleValidation.aggregate_array('predicted'));
// Compute residuals
var residualsValidation = observationValidation.subtract(predictionValidation);
// Compute RMSE with equation, print it
var rmseValidation = residualsValidation.pow(2).reduce('mean', [0]).sqrt();
print('Validation RMSE', rmseValidation);

圖片

導出

// **** Export regression **** //
// create file name for export
var exportName = 'DSM_RF_regression';// If using gmail: Export to Drive
Export.image.toDrive({image: regression,description: exportName,folder: "DigitalSoilMapping",fileNamePrefix: exportName,region: roi,scale: 30,crs: crs,maxPixels: 1e13
});

運行腳本后,窗格右側的“任務”選項卡將變為橙色,表示有可以運行的導出任務。

圖片

單擊任務選項卡。

在窗格中找到相應的任務,然后單擊“運行”。

圖片

在彈出的窗口中,您將看到導出參數。我們已經在代碼中指定了這些。請再次檢查以確保一切正常,然后單擊“運行”。導出可能需要 10 分鐘以上才能完成,所以請耐心等待!

圖片

請注意,我們正在將導出文件發送到您的 Google Drive 中名為“DigitalSoilMapping”的文件夾中。如果此文件夾尚不存在,則此命令將創建它。

導航到您的 google 驅動器,找到 DigitalSoilMapping 文件夾,然后單擊將其打開。

右鍵單擊下載文件,該文件的標題應為“Essex_VT_DSM_regression.tif”。現在,您可以在您選擇的 GIS 中打開它。

恭喜!您已成功完成此練習。

// 以下是本次實驗的全部代碼
var VT_boundary = ee.FeatureCollection("projects/ee-yelu/assets/VT_boundary"),VT_pedons = ee.FeatureCollection("projects/ee-yelu/assets/essex_pedons_all");
// 研究區
var roi = VT_boundary.union(); 
Map.addLayer(roi, {}, 'shp', false);
var crs = 'EPSG:4326'; // EPSG number for output projection. 32618 = WGS84/UTM Zone 18N. For more info- http://spatialreference.org/ref/epsg/
// S2 影像去云與合成
function maskS2clouds(image) {var qa = image.select('QA60');// Bits 10 and 11 are clouds and cirrus, respectively.var cloudBitMask = 1 << 10;var cirrusBitMask = 1 << 11;// Both flags should be set to zero, indicating clear conditions.var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));return image.updateMask(mask).divide(10000);
}// 計算特征
// NDVI
function NDVI(img) {var ndvi = img.expression("(NIR-R)/(NIR+R)", {"R": img.select("B4"),"NIR": img.select("B8"),});return img.addBands(ndvi.rename("NDVI"));
}// NDWI
function NDWI(img) {var ndwi = img.expression("(G-MIR)/(G+MIR)", {"G": img.select("B3"),"MIR": img.select("B8"),});return img.addBands(ndwi.rename("NDWI"));
}// NDBI
function NDBI(img) {var ndbi = img.expression("(SWIR-NIR)/(SWIR-NIR)", {"NIR": img.select("B8"),"SWIR": img.select("B12"),});return img.addBands(ndbi.rename("NDBI"));
}//SAVI
function SAVI(image) {var savi = image.expression('(NIR - RED) * (1 + 0.5)/(NIR + RED + 0.5)', {'NIR': image.select('B8'),'RED': image.select('B4')}).float();return image.addBands(savi.rename('SAVI'));
}//IBI 
function IBI(image) {// Add Index-Based Built-Up Index (IBI)var ibiA = image.expression('2 * SWIR1 / (SWIR1 + NIR)', {'SWIR1': image.select('B6'),'NIR': image.select('B5')}).rename(['IBI_A']);var ibiB = image.expression('(NIR / (NIR + RED)) + (GREEN / (GREEN + SWIR1))', {'NIR': image.select('B8'),'RED': image.select('B4'),'GREEN': image.select('B3'),'SWIR1': image.select('B11')}).rename(['IBI_B']);var ibiAB = ibiA.addBands(ibiB);var ibi = ibiAB.normalizedDifference(['IBI_A', 'IBI_B']);return image.addBands(ibi.rename(['IBI']));
}//RVI
function RVI(image) {var rvi = image.expression('NIR/Red', {'NIR': image.select('B8'),'Red': image.select('B4')});return image.addBands(rvi.rename('RVI'));
}//DVI
function DVI(image) {var dvi = image.expression('NIR - Red', {'NIR': image.select('B8'),'Red': image.select('B4')}).float();return image.addBands(dvi.rename('DVI'));
}// 創建var inStack = ee.Image()
for (var i = 0; i < 3; i += 1) {var start = ee.Date('2019-03-01').advance(30 * i, 'day');print(start)var end = start.advance(30, 'day');var dataset = ee.ImageCollection('COPERNICUS/S2_SR').filterDate(start, end).filterBounds(roi).filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',75)).map(maskS2clouds).map(NDVI).map(NDWI).map(NDBI).map(SAVI).map(IBI).map(RVI).map(DVI);var inStack_monthly = dataset.median().clip(roi);// 可視化var visualization = {min: 0.0,max: 0.3,bands: ['B4', 'B3', 'B2'],};Map.addLayer(inStack_monthly, visualization, 'S2_' + i, false);// // 紋理特征var B8 = inStack_monthly.select('B8').multiply(100).toInt16();var glcm = B8.glcmTexture({size: 3});var contrast = glcm.select('B8_contrast');var var_ = glcm.select('B8_var');var savg = glcm.select('B8_savg');var dvar = glcm.select('B8_dvar');inStack_monthly = inStack_monthly.addBands([contrast, var_, savg, dvar])// Sentinel-1 var imgVV = ee.ImageCollection('COPERNICUS/S1_GRD').filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV')).filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH')).filter(ee.Filter.eq('instrumentMode', 'IW')).filterBounds(roi).map(function(image) {var edge = image.lt(-30.0);var maskedImage = image.mask().and(edge.not());return image.updateMask(maskedImage);});var img_S1_asc = imgVV.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));var VV_img = img_S1_asc.filterDate(start, end).select("VV").median().clip(roi);var VH_img = img_S1_asc.filterDate(start, end).select("VH").median().clip(roi);inStack_monthly = inStack_monthly.addBands([VV_img, VH_img])print("img_S2_monthly", inStack_monthly)inStack = inStack.addBands(inStack_monthly.select(inStack_monthly.bandNames()))
}
inStack = inStack.select(inStack.bandNames().remove("constant"))
print(inStack)
inStack.reproject(crs, null, 30);
Map.centerObject(inStack)var trainingFeatureCollection = ee.FeatureCollection(VT_pedons, 'geometry');
print(bands)
var training = inStack.reduceRegions({collection:trainingFeatureCollection,reducer :ee.Reducer.mean(),scale:100,tileScale :5,crs: crs
});
print(training)
var bands = inStack.bandNames()var trainging2 = training.filter(ee.Filter.notNull(bands));
// print(trainging2);// 在training要素集中增加一個random屬性,值為0到1的隨機數
var withRandom = trainging2.randomColumn({columnName:'random',seed:2,distribution: 'uniform',
});var split = 0.8; 
var trainingData = withRandom.filter(ee.Filter.lt('random', split));
var validationData = withRandom.filter(ee.Filter.gte('random', split));    // train the RF classification model
var classifier = ee.Classifier.smileRandomForest(100, null, 1, 0.5, null, 0).setOutputMode('REGRESSION').train({features: trainingData, classProperty: 'REDOX_CM', inputProperties: bands});
print(classifier);// Classify the input imagery.
var regression = inStack.classify(classifier, 'predicted');
print(regression);var palettes = require('users/gena/packages:palettes');
var palette = palettes.crameri.nuuk[25];
print("regression", regression)
// Display the input imagery and the regression classification.
// get dictionaries of min & max predicted value
var regressionMin = (regression.reduceRegion({reducer: ee.Reducer.min(),scale: 100,crs: crs,bestEffort: true,tileScale: 5,geometry: roi
}));
var regressionMax = (regression.reduceRegion({reducer: ee.Reducer.max(),scale: 100,crs: crs,bestEffort: true,tileScale: 5,geometry: roi
}));// Add to map
var viz = {palette: palette,min: regressionMin.getNumber('predicted').getInfo(),max: regressionMax.getNumber('predicted').getInfo()
};
Map.addLayer(regression, viz, 'Regression');// Create the panel for the legend items.
var legend = ui.Panel({style: {position: 'bottom-left',padding: '8px 15px'}
});// Create and add the legend title.
var legendTitle = ui.Label({value: 'Legend',style: {fontWeight: 'bold',fontSize: '18px',margin: '0 0 4px 0',padding: '0'}
});
legend.add(legendTitle);// create the legend image
var lon = ee.Image.pixelLonLat().select('latitude');
var gradient = lon.multiply((viz.max - viz.min) / 100.0).add(viz.min);
var legendImage = gradient.visualize(viz);// create text on top of legend
var panel = ui.Panel({widgets: [ui.Label(viz['max'])],
});legend.add(panel);// create thumbnail from the image
var thumbnail = ui.Thumbnail({image: legendImage,params: {bbox: '0,0,10,100',dimensions: '10x200'},style: {padding: '1px',position: 'bottom-center'}
});// add the thumbnail to the legend
legend.add(thumbnail);// create text on top of legend
var panel = ui.Panel({widgets: [ui.Label(viz['min'])],
});legend.add(panel);
Map.add(legend);// Zoom to the regression on the map
Map.centerObject(roi, 11);// Get variable importance
var dict = classifier.explain();
print("Classifier information:", dict);
var variableImportance = ee.Feature(null, ee.Dictionary(dict).get('importance'));
// Make chart, print it
var chart =ui.Chart.feature.byProperty(variableImportance).setChartType('ColumnChart').setOptions({title: 'Random Forest Variable Importance',legend: {position: 'none'},hAxis: {title: 'Bands'},vAxis: {title: 'Importance'}});
print(chart);// Set chart options
var options = {lineWidth: 1,pointSize: 2,hAxis: {title: 'Redox (cm)'},vAxis: {title: 'Num of pixels'},title: 'Number of pixels by redox depth'
};
var regressionPixelChart = ui.Chart.image.histogram({image: ee.Image(regression),region: roi,scale:100}).setOptions(options);
print(regressionPixelChart);// Get predicted regression points in same location as training data
var predictedTraining = regression.sampleRegions({collection: trainingData,scale:100,geometries: true,projection:crs
});
// Separate the observed (REDOX_CM) and predicted (regression) properties
var sampleTraining = predictedTraining.select(['REDOX_CM', 'predicted']);
// Create chart, print it
var chartTraining = ui.Chart.feature.byFeature({features:sampleTraining, xProperty:'REDOX_CM', yProperties:['predicted']}).setChartType('ScatterChart').setOptions({title: 'Predicted vs Observed - Training data ',hAxis: {'title': 'observed'},vAxis: {'title': 'predicted'},pointSize: 3,trendlines: {0: {showR2: true,visibleInLegend: true},1: {showR2: true,visibleInLegend: true}}});
print(chartTraining);// **** Compute RSME **** 
// Get array of observation and prediction values 
var observationTraining = ee.Array(sampleTraining.aggregate_array('REDOX_CM'));
var predictionTraining = ee.Array(sampleTraining.aggregate_array('predicted'));
print('observationTraining', observationTraining)
print('predictionTraining', predictionTraining)// Compute residuals
var residualsTraining = observationTraining.subtract(predictionTraining);
print('residualsTraining', residualsTraining)
// Compute RMSE with equation, print it
var rmseTraining = residualsTraining.pow(2).reduce('mean', [0]).sqrt();
print('Training RMSE', rmseTraining);// Get predicted regression points in same location as validation data
var predictedValidation = regression.sampleRegions({collection: validationData,scale:100,geometries: true
});
// Separate the observed (REDOX_CM) and predicted (regression) properties
var sampleValidation = predictedValidation.select(['REDOX_CM', 'predicted']);
// Create chart, print it
var chartValidation = ui.Chart.feature.byFeature(sampleValidation, 'predicted', 'REDOX_CM').setChartType('ScatterChart').setOptions({title: 'Predicted vs Observed - Validation data',hAxis: {'title': 'predicted'},vAxis: {'title': 'observed'},pointSize: 3,trendlines: {0: {showR2: true,visibleInLegend: true},1: {showR2: true,visibleInLegend: true}}});
print(chartValidation);// Get array of observation and prediction values 
var observationValidation = ee.Array(sampleValidation.aggregate_array('REDOX_CM'));
var predictionValidation = ee.Array(sampleValidation.aggregate_array('predicted'));
// Compute residuals
var residualsValidation = observationValidation.subtract(predictionValidation);
// Compute RMSE with equation, print it
var rmseValidation = residualsValidation.pow(2).reduce('mean', [0]).sqrt();
print('Validation RMSE', rmseValidation);// **** Export regression **** //
// create file name for export
var exportName = 'DSM_RF_regression';// If using gmail: Export to Drive
Export.image.toDrive({image: regression,description: exportName,folder: "DigitalSoilMapping",fileNamePrefix: exportName,region: roi,scale: 30,crs: crs,maxPixels: 1e13
});

本文來自互聯網用戶投稿,該文觀點僅代表作者本人,不代表本站立場。本站僅提供信息存儲空間服務,不擁有所有權,不承擔相關法律責任。
如若轉載,請注明出處:http://www.pswp.cn/news/207307.shtml
繁體地址,請注明出處:http://hk.pswp.cn/news/207307.shtml
英文地址,請注明出處:http://en.pswp.cn/news/207307.shtml

如若內容造成侵權/違法違規/事實不符,請聯系多彩編程網進行投訴反饋email:809451989@qq.com,一經查實,立即刪除!

相關文章

Codeforces Round 851 (Div. 2 D:枚舉+組合 Edp)

A - One and Two 相當于找第一個位置前后2的個數相同 #include<bits/stdc.h> using namespace std; const int N 1e610,mod998244353; #define int long long typedef long long LL; typedef pair<int, int> PII; const long long inf1e17; int n,m,k; int a[N]…

有哪些值得分享的銷售拓客技巧?

拓客對于銷售的重要性 拓客&#xff08;Toker&#xff09;是一個商業上的名詞&#xff0c;核心就是提高售前服務、市場推廣的水平&#xff0c;從而挖掘出潛在客戶的隱形需求&#xff08;或稱軟需求&#xff09;。 拓客的核心&#xff0c;其實就是提高售前服務、市場推廣的水平…

如何部署自己的服務渲染頁面為Pdf文檔

前言 相信大家都覺得官方發布的文檔生成模塊https://docs.mendix.com/appstore/modules/document-generation/很有用&#xff0c;它能把Mendix頁面像素級導出到Pdf文件中&#xff0c;這對于歸檔等業務非常有價值。但部署依賴公有云提供的渲染服務&#xff0c;而中國本土用戶對…

折半查找(數據結構實訓)

題目&#xff1a; 標準輸入輸出 題目描述&#xff1a; 實現折半查找。要求查找給定的值在數據表中相應的存儲位置。本題目假定輸入元素均按非降序輸入。 輸入&#xff1a; 輸入包含若干個測試用例&#xff0c;第一行為測試用例個數k。每個測試用例占3行&#xff0c;其中第一行為…

初識人工智能,一文讀懂過擬合欠擬合和模型壓縮的知識文集(3)

&#x1f3c6;作者簡介&#xff0c;普修羅雙戰士&#xff0c;一直追求不斷學習和成長&#xff0c;在技術的道路上持續探索和實踐。 &#x1f3c6;多年互聯網行業從業經驗&#xff0c;歷任核心研發工程師&#xff0c;項目技術負責人。 &#x1f389;歡迎 &#x1f44d;點贊?評論…

SQL存儲過程和視圖

1 存儲過程 存儲過程是事先編寫好、存儲在數據庫中的一組SQL命令集合。用來完成對數據庫的指定操作。 1.1 優缺點 優點&#xff1a; 1&#xff09;提高系統性能。創建時進行編譯&#xff0c;隨后存放在數據庫服務器的過程高速緩存中&#xff0c;之后不需要再次執行分析和編…

uniapp app將base64保存到相冊,uniapp app將文件流保存到相冊

如果是文件流可以先轉base64詳情見>uniapp 顯示文件流圖片-CSDN博客 onDown(){let base64 this.qrcodeUrl ; // base64地址const bitmap new plus.nativeObj.Bitmap("test");bitmap.loadBase64Data(base64, function() {const url "_doc/" new Dat…

Backend - Dbeaver

目錄 一、說明 二、下載并安裝 &#xff08;一&#xff09;官網下載 &#xff08;二&#xff09;安裝 三、使用 &#xff08;一&#xff09;操作步驟 &#xff08;二&#xff09;相關問題&#xff1a;無法加載驅動類oracle.jdbc.oracledriver 1. 新建驅動 2. 再重新連接數據庫 …

PyTorch2.0環境搭建

一、安裝python并配置環境變量 1、打開python官網&#xff0c;下載并安裝 Welcome to Python.org 下載 尋找版本&#xff1a;推薦使用3.9版本&#xff0c;或其他表中顯示為安全&#xff08;security&#xff09;的版本 安裝&#xff1a;&#xff08;略&#xff09; 2、配置環…

數據增強改進,實現檢測目標copypaste,增加目標數據量,提升精度

???YOLOv8實戰寶典--星級指南:從入門到精通,您不可錯過的技巧 ??-- 聚焦于YOLO的 最新版本, 對頸部網絡改進、添加局部注意力、增加檢測頭部,實測漲點 ?? 深入淺出YOLOv8:我的專業筆記與技術總結 ??-- YOLOv8輕松上手, 適用技術小白,文章代碼齊全,僅需 …

python圣誕樹代碼編程

以下是一個簡單的Python圣誕樹代碼&#xff1a; def draw_tree(height): for i in range(height): print( * (height - i - 1) * * (2 * i 1)) print( * (height - 1) |)draw_tree(10) 這個函數會繪制一個等腰三角形&#xff0c;其中每一行的星號數量從1開…

Java基礎知識

JVM&#xff0c;JRE&#xff0c;JDK JVM 運行Java字節碼的機器 JRE Java運行時環境&#xff0c;包括JVM&#xff0c;Java類庫&#xff0c;運行時類庫&#xff0c;國際化支持&#xff0c;安全管理器&#xff0c;啟動器等 比JVM多的內容 Java類庫&#xff1a;提供大量已經實…

【TiDB理論知識09】TiFlash

一 TiFlash架構 二 TiFlash 核心特性 TiFlash 主要有 異步復制、一致性、智能選擇、計算加速 等幾個核心特性。 1 異步復制 TiFlash 中的副本以特殊角色 (Raft Learner) 進行異步的數據復制&#xff0c;這表示當 TiFlash 節點宕機或者網絡高延遲等狀況發生時&#xff0c;Ti…

億勝盈科ATR2037 無限射頻前端低噪聲放大器

億勝盈科ATR2037 是一款應用于無線通信射頻前端&#xff0c;工作頻段為 0.7 到 6GHz 的超低噪聲放大器。 ATR2037 低噪聲放大器采用先進的 GaAs pHEMT 工藝設計和制作&#xff0c;ATR2037 低噪聲放大器在整個工作頻段內可以獲得非常好的射頻性能超低噪聲系數。 億勝盈科ATR203…

WGCLOUD v3.5.0 新增支持監測交換機的接口狀態UP DOWN

WGCLOUD v3.5.0開始 可以監測交換機或SNMP設備的接口狀態了&#xff0c;直接上圖

什么是ElasticSearch中的過濾器?

在Elasticsearch中&#xff0c;過濾器&#xff08;Filters&#xff09;是一種用于在查詢中篩選文檔的強大工具。過濾器可以根據特定條件來評估文檔是否符合搜索查詢。這些條件通常應用于字段數據&#xff0c;并根據匹配結果返回符合條件的文檔。 過濾器的主要優點包括&#xf…

如何給網頁和代碼做HTML加密?

? 本篇文章給大家談談html混淆加密在線&#xff0c;以及HTML在線加密對應的知識點&#xff0c;希望對各位有所幫助&#xff0c;不要忘了收藏本站喔。 如何給代碼加密? 1、源代碼加密軟件推薦使用德人合科技的加密軟件&#xff0c;是一套從源頭上保障數據安全和使用安全的軟…

vue2+datav可視化數據大屏(1)

開始 &#x1f4d3; 最近打算出一個前端可視化數據大屏的系列專欄&#xff0c;這次將很全面的教大家設計可視化大屏&#xff0c;從開始到打包結束&#xff0c;其中&#xff0c;包括如何設計框架&#xff0c;如何封裝axios&#xff0c;等等&#xff0c;本次使用的數據均為mock數…

linux C++監聽管道文件方式

方式一&#xff08;傳統讀取文件&#xff0c;一直監聽循環讀取文件&#xff09; 非阻塞打開文件&#xff0c;用read循環定時讀取&#xff0c;性能不好 代碼如下&#xff1a; #include <iostream> #include <fstream> #include <functional> #include <…

spring boot項目如何自定義參數校驗規則

spring boot項目對參數進行校驗時&#xff0c;比如非空校驗&#xff0c;可以直接用validation包里面自帶的注解。但是對于一些復雜的參數校驗&#xff0c;自帶的校驗規則無法滿足要求&#xff0c;此時需要我們自定義參數校驗規則。自定義校驗規則和自帶的規則實現方式一樣&…