ฟิลด์การฝังสามารถใช้เป็นอินพุต/ตัวทำนายฟีเจอร์สำหรับการถดถอยในลักษณะเดียวกับที่ใช้สำหรับการจัดประเภท
ในบทแนะนำนี้ เราจะมาดูวิธีใช้เลเยอร์ฟิลด์การฝัง 64D เป็นอินพุตในการวิเคราะห์การถดถอยเชิงพหุเพื่อคาดการณ์ชีวมวลเหนือพื้นดิน (AGB)
ภารกิจ Global Ecosystem Dynamics Investigation (GEDI) ของ NASA รวบรวมการวัด LIDAR ตามแนวตัดขวางภาคพื้นดินที่ความละเอียดเชิงพื้นที่ 30 ม. ที่ช่วง 60 ม. เราจะใช้ชุดข้อมูล GEDI L4A Raster Aboveground Biomass Density ซึ่งมีค่าประมาณแบบจุดของความหนาแน่นของมวลชีวภาพเหนือพื้นดิน (AGBD) ที่จะใช้เป็นตัวแปรที่คาดการณ์ไว้ในโมเดลการถดถอย
เลือกภูมิภาค
มาเริ่มด้วยการกำหนดภูมิภาคที่สนใจกัน สำหรับบทแนะนำนี้ เราจะเลือกภูมิภาคในเทือกเขา Western Ghats ของอินเดีย และกำหนดรูปหลายเหลี่ยมเป็นตัวแปรเรขาคณิต หรือจะใช้เครื่องมือวาดในตัวแก้ไขโค้ดเพื่อวาดรูปหลายเหลี่ยมรอบๆ ภูมิภาคที่สนใจซึ่งจะบันทึกเป็นตัวแปรเรขาคณิตในการนำเข้าก็ได้ นอกจากนี้ เรายังใช้แผนที่ฐานจากดาวเทียม ซึ่งช่วยให้ค้นหาพื้นที่ที่มีพืชพรรณได้ง่าย
var geometry = ee.Geometry.Polygon([[
[74.322, 14.981],
[74.322, 14.765],
[74.648, 14.765],
[74.648, 14.980]
]]);
// Use the satellite basemap
Map.setOptions('SATELLITE');
รูปภาพ: การเลือกพื้นที่ที่สนใจสำหรับการคาดการณ์มวลชีวภาพเหนือพื้นดิน
เลือกระยะเวลา
เลือกปีที่ต้องการเรียกใช้การถดถอย โปรดทราบว่าการฝังดาวเทียมจะรวบรวมข้อมูลเป็นช่วงปี ดังนั้นเราจึงกำหนดระยะเวลาตลอดทั้งปี
var startDate = ee.Date.fromYMD(2022, 1, 1);
var endDate = startDate.advance(1, 'year');
เตรียมชุดข้อมูลการฝังดาวเทียม
ระบบจะใช้รูปภาพการฝังดาวเทียมแบบ 64 แบนด์เป็นตัวทำนายสำหรับการถดถอย เราโหลดชุดข้อมูลการฝังดาวเทียม กรองรูปภาพสำหรับปีและภูมิภาคที่เลือก
var embeddings = ee.ImageCollection('GOOGLE/SATELLITE_EMBEDDING/V1/ANNUAL');
var embeddingsFiltered = embeddings
.filter(ee.Filter.date(startDate, endDate))
.filter(ee.Filter.bounds(geometry));
รูปภาพการฝังดาวเทียมจะแบ่งเป็นตารางในไทล์และแสดงในการฉายภาพสำหรับโซน UTM ของไทล์ ด้วยเหตุนี้ เราจึงได้ไทล์การฝังสภาพพื้นผิวจากดาวเทียมหลายรายการที่ครอบคลุมภูมิภาคที่สนใจ หากต้องการได้ภาพเดียว เราต้องต่อภาพเหล่านั้น ใน Earth Engine ระบบจะกำหนดการฉายภาพเริ่มต้นให้กับภาพอินพุตที่ต่อกันเป็นผืน ซึ่งก็คือ WGS84 ที่มีมาตราส่วน 1 องศา เนื่องจากเราจะรวบรวมและฉายภาพโมเสกนี้อีกครั้งในภายหลังในบทแนะนำ การคงการฉายภาพต้นฉบับไว้จึงเป็นประโยชน์ เราสามารถดึงข้อมูลการฉายจากไทล์รายการใดรายการหนึ่งและตั้งค่าในโมเสกได้โดยใช้ฟังก์ชัน setDefaultProjection()
// Extract the projection of the first band of the first image
var embeddingsProjection = ee.Image(embeddingsFiltered.first()).select(0).projection();
// Set the projection of the mosaic to the extracted projection
var embeddingsImage = embeddingsFiltered.mosaic()
.setDefaultProjection(embeddingsProjection);
เตรียมภาพโมเสก GEDI L4A
เนื่องจากระบบจะใช้ค่าประมาณชีวมวล GEDI เพื่อฝึกโมเดลการถดถอยของเรา จึงจำเป็นอย่างยิ่งที่จะต้องกรองข้อมูล GEDI ที่ไม่ถูกต้องหรือไม่น่าเชื่อถือก่อนนำไปใช้ เราใช้มาสก์หลายรายการเพื่อนำการวัดที่อาจมีข้อผิดพลาดออก
- นำการวัดทั้งหมดที่ไม่เป็นไปตามข้อกำหนดด้านคุณภาพออก (l4_quality_flag = 0 และ degrade_flag > 0)
- นำการวัดทั้งหมดที่มีข้อผิดพลาดสัมพัทธ์สูง ("agbd_se" / "agbd" > 50%) ออก
- นำการวัดทั้งหมดบนความชัน > 30% ออกโดยอิงตามรูปแบบความสูงเชิงเลข (DEM) ของ Copernicus GLO-30
สุดท้าย เราจะเลือกการวัดที่เหลือทั้งหมดสำหรับช่วงเวลาและภูมิภาคที่สนใจ แล้วสร้างภาพโมเสก
var gedi = ee.ImageCollection('LARSE/GEDI/GEDI04_A_002_MONTHLY');
// Function to select the highest quality GEDI data
var qualityMask = function(image) {
return image.updateMask(image.select('l4_quality_flag').eq(1))
.updateMask(image.select('degrade_flag').eq(0));
};
// Function to mask unreliable GEDI measurements
// with a relative standard error > 50%
// agbd_se / agbd > 0.5
var errorMask = function(image) {
var relative_se = image.select('agbd_se')
.divide(image.select('agbd'));
return image.updateMask(relative_se.lte(0.5));
};
// Function to mask GEDI measurements on slopes > 30%
var slopeMask = function(image) {
// Use Copernicus GLO-30 DEM for calculating slope
var glo30 = ee.ImageCollection('COPERNICUS/DEM/GLO30');
var glo30Filtered = glo30
.filter(ee.Filter.bounds(geometry))
.select('DEM');
// Extract the projection
var demProj = glo30Filtered.first().select(0).projection();
// The dataset consists of individual images
// Create a mosaic and set the projection
var elevation = glo30Filtered.mosaic().rename('dem')
.setDefaultProjection(demProj);
// Compute the slope
var slope = ee.Terrain.slope(elevation);
return image.updateMask(slope.lt(30));
};
var gediFiltered = gedi
.filter(ee.Filter.date(startDate, endDate))
.filter(ee.Filter.bounds(geometry));
var gediProjection = ee.Image(gediFiltered.first())
.select('agbd').projection();
var gediProcessed = gediFiltered
.map(qualityMask)
.map(errorMask)
.map(slopeMask);
var gediMosaic = gediProcessed.mosaic()
.select('agbd').setDefaultProjection(gediProjection);
// Visualize the GEDI Mosaic
var gediVis = {
min: 0,
max: 200,
palette: ['#edf8fb','#b2e2e2','#66c2a4','#2ca25f','#006d2c'],
bands: ['agbd']
};
Map.addLayer(gediMosaic, gediVis, 'GEDI L4A (Filtered)', false);
สุ่มตัวอย่างใหม่และรวบรวมอินพุต
ก่อนที่จะสุ่มตัวอย่างพิกเซลเพื่อฝึกโมเดลการถดถอย เราจะสุ่มตัวอย่างใหม่และฉายอินพุตใหม่ไปยังตารางพิกเซลเดียวกัน การวัดของ GEDI มีความแม่นยำในแนวนอน +/- 9 ม. ซึ่งจะเกิดปัญหาเมื่อจับคู่ค่า GEDI AGB กับพิกเซลการฝังดาวเทียม เราจึงทำการสุ่มตัวอย่างใหม่และรวบรวมรูปภาพอินพุตทั้งหมดลงในตารางกริดพิกเซลที่ใหญ่ขึ้นโดยใช้ค่าเฉลี่ยจากพิกเซลต้นฉบับ ซึ่งยังช่วยลดสัญญาณรบกวนจากข้อมูลและช่วยสร้างโมเดลแมชชีนเลิร์นนิงที่ดีขึ้นด้วย
// Choose the grid size and projection
var gridScale = 100;
var gridProjection = ee.Projection('EPSG:3857')
.atScale(gridScale);
// Create a stacked image with predictor and predicted variables
var stacked = embeddingsImage.addBands(gediMosaic);
// Set the resampling mode
var stacked = stacked.resample('bilinear');
// Aggregate pixels with 'mean' statistics
var stackedResampled = stacked
.reduceResolution({
reducer: ee.Reducer.mean(),
maxPixels: 1024
})
.reproject({
crs: gridProjection
});
// As larger GEDI pixels contain masked original
// pixels, it has a transparency mask.
// We update the mask to remove the transparency
var stackedResampled = stackedResampled
.updateMask(stackedResampled.mask().gt(0));
การฉายซ้ำและการรวมพิกเซลเป็นการดำเนินการที่มีค่าใช้จ่ายสูง และแนวทางปฏิบัติแนะนำคือการส่งออกรูปภาพที่ซ้อนกันเป็นชิ้นงาน แล้วใช้รูปภาพที่คำนวณไว้ล่วงหน้าในขั้นตอนต่อๆ ไป ซึ่งจะช่วยแก้ไขข้อผิดพลาดการคำนวณหมดเวลาหรือหน่วยความจำของผู้ใช้เกินเมื่อทำงานกับภูมิภาคขนาดใหญ่
// Replace this with your asset folder
// The folder must exist before exporting
var exportFolder = 'projects/spatialthoughts/assets/satellite_embedding/';
var mosaicExportImage = 'gedi_mosaic';
var mosaicExportImagePath = exportFolder + mosaicExportImage;
Export.image.toAsset({
image: stackedResampled.clip(geometry),
description: 'GEDI_Mosaic_Export',
assetId: mosaicExportImagePath,
region: geometry,
scale: gridScale,
maxPixels: 1e10
});
เริ่มงานส่งออกและรอจนกว่าจะเสร็จสิ้น เมื่อเสร็จแล้ว เราจะนำเข้าชิ้นงานและสร้างโมเดลต่อไป
// Use the exported asset
var stackedResampled = ee.Image(mosaicExportImagePath);
ดึงข้อมูลฟีเจอร์การฝึก
เรามีข้อมูลอินพุตพร้อมสำหรับการแยกฟีเจอร์การฝึกแล้ว เราใช้แถบการฝังดาวเทียมเป็นตัวแปรตาม (ตัวทำนาย) และค่า AGBD ของ GEDI เป็นตัวแปรอิสระ (ค่าที่คาดการณ์) ในโมเดลการถดถอย เราสามารถดึงค่าที่ตรงกันในแต่ละพิกเซลและเตรียมชุดข้อมูลการฝึกได้ รูปภาพ GEDI ของเราส่วนใหญ่จะถูกมาสก์และมีค่าในพิกเซลเพียงส่วนย่อยๆ เท่านั้น หากเราใช้ sample() ระบบจะแสดงค่าว่างเป็นส่วนใหญ่ เราจึงสร้างแถบชั้นเรียนจากมาสก์ GEDI และใช้ stratifiedSample() เพื่อให้มั่นใจว่าเราจะสุ่มตัวอย่างจากพิกเซลที่ไม่ได้มาสก์
var predictors = embeddingsImage.bandNames();
var predicted = gediMosaic.bandNames().get(0);
print('predictors', predictors);
print('predicted', predicted);
var predictorImage = stackedResampled.select(predictors);
var predictedImage = stackedResampled.select([predicted]);
var classMask = predictedImage.mask().toInt().rename('class');
var numSamples = 1000;
// We set classPoints to [0, numSamples]
// This will give us 0 points for class 0 (masked areas)
// and numSample points for class 1 (non-masked areas)
var training = stackedResampled.addBands(classMask)
.stratifiedSample({
numPoints: numSamples,
classBand: 'class',
region: geometry,
scale: gridScale,
classValues: [0, 1],
classPoints: [0, numSamples],
dropNulls: true,
tileScale: 16,
});
print('Number of Features Extracted', training.size());
print('Sample Training Feature', training.first());
ฝึกโมเดลการถดถอย
ตอนนี้เราพร้อมที่จะฝึกโมเดลแล้ว ตัวแยกประเภทจำนวนมากใน Earth Engine สามารถใช้สำหรับงานทั้งการแยกประเภทและการถดถอย เนื่องจากเราต้องการคาดการณ์ค่าตัวเลข (แทนที่จะเป็นคลาส) เราจึงตั้งค่าเครื่องมือแยกประเภทให้ทำงานในโหมด REGRESSION และฝึกโดยใช้ข้อมูลการฝึกได้ เมื่อฝึกโมเดลแล้ว เราจะเปรียบเทียบการคาดการณ์ของโมเดลกับค่าอินพุต และคำนวณค่าความคลาดเคลื่อนกำลังสองเฉลี่ย (rmse) และสัมประสิทธิ์สหสัมพันธ์ r^2 เพื่อตรวจสอบประสิทธิภาพของโมเดลได้
// Use the RandomForest classifier and set the
// output mode to REGRESSION
var model = ee.Classifier.smileRandomForest(50)
.setOutputMode('REGRESSION')
.train({
features: training,
classProperty: predicted,
inputProperties: predictors
});
// Get model's predictions for training samples
var predicted = training.classify({
classifier: model,
outputName: 'agbd_predicted'
});
// Calculate RMSE
var calculateRmse = function(input) {
var observed = ee.Array(
input.aggregate_array('agbd'));
var predicted = ee.Array(
input.aggregate_array('agbd_predicted'));
var rmse = observed.subtract(predicted).pow(2)
.reduce('mean', [0]).sqrt().get([0]);
return rmse;
};
var rmse = calculateRmse(predicted);
print('RMSE', rmse);
// Create a plot of observed vs. predicted values
var chart = ui.Chart.feature.byFeature({
features: predicted.select(['agbd', 'agbd_predicted']),
xProperty: 'agbd',
yProperties: ['agbd_predicted'],
}).setChartType('ScatterChart')
.setOptions({
title: 'Aboveground Biomass Density (Mg/Ha)',
dataOpacity: 0.8,
hAxis: {'title': 'Observed'},
vAxis: {'title': 'Predicted'},
legend: {position: 'right'},
series: {
0: {
visibleInLegend: false,
color: '#525252',
pointSize: 3,
pointShape: 'triangle',
},
},
trendlines: {
0: {
type: 'linear',
color: 'black',
lineWidth: 1,
pointSize: 0,
labelInLegend: 'Linear Fit',
visibleInLegend: true,
showR2: true
}
},
chartArea: {left: 100, bottom: 100, width: '50%'},
});
print(chart);
รูปภาพ: ค่า AGBD ที่สังเกตได้เทียบกับค่าที่โมเดลคาดการณ์
สร้างการคาดการณ์สำหรับค่าที่ไม่รู้จัก
เมื่อพอใจกับโมเดลแล้ว เราจะใช้โมเดลที่ฝึกแล้วเพื่อสร้างการคาดการณ์ในตำแหน่งที่ไม่รู้จักจากรูปภาพที่มีแถบตัวแปรทำนาย
// We set the band name of the output image as 'agbd'
var predictedImage = stackedResampled.classify({
classifier: model,
outputName: 'agbd'
});
ตอนนี้รูปภาพที่มีค่า AGBD ที่คาดการณ์ไว้ในแต่ละพิกเซลพร้อมสำหรับการส่งออกแล้ว เราจะใช้ข้อมูลนี้ในส่วนถัดไปเพื่อแสดงผลลัพธ์เป็นภาพ
// Replace this with your asset folder
// The folder must exist before exporting
var exportFolder = 'projects/spatialthoughts/assets/satellite_embedding/';
var predictedExportImage = 'predicted_agbd';
var predictedExportImagePath = exportFolder + predictedExportImage;
Export.image.toAsset({
image: predictedImage.clip(geometry),
description: 'Predicted_Image_Export',
assetId: predictedExportImagePath,
region: geometry,
scale: gridScale,
maxPixels: 1e10
});
เริ่มงานส่งออกและรอจนกว่าจะเสร็จสิ้น เมื่อเสร็จแล้ว เราจะนำเข้าชิ้นงานและแสดงผลลัพธ์
var predictedImage = ee.Image(predictedExportImagePath);
// Visualize the image
var gediVis = {
min: 0,
max: 200,
palette: ['#edf8fb','#b2e2e2','#66c2a4','#2ca25f','#006d2c'],
bands: ['agbd']
};
Map.addLayer(predictedImage, gediVis, 'Predicted AGBD');
รูปภาพ: AGBD ที่คาดการณ์ สีเขียวเข้มขึ้นแสดงถึงความหนาแน่นของมวลชีวภาพที่คาดการณ์ไว้มากขึ้น
ประมาณมวลชีวภาพทั้งหมด
ตอนนี้เรามีค่า AGBD ที่คาดการณ์ไว้สำหรับแต่ละพิกเซลของรูปภาพแล้ว และสามารถใช้เพื่อประมาณปริมาณชีวมวลทั้งหมดเหนือพื้นดิน (AGB) ในภูมิภาคได้ แต่เราต้องนำพิกเซลทั้งหมดที่เป็นของพื้นที่ที่ไม่มีพืชออกก่อน เราสามารถใช้ชุดข้อมูลการปกคลุมของพื้นดิน ESA WorldCover และเลือกพิกเซลที่มีพืช
// GEDI data is processed only for certain landcovers
// from Plant Functional Types (PFT) classification
// https://doi.org/10.1029/2022EA002516
// Here we use ESA WorldCover v200 product to
// select landcovers representing vegetated areas
var worldcover = ee.ImageCollection('ESA/WorldCover/v200').first();
// Aggregate pixels to the same grid as other dataset
// with 'mode' value.
// i.e. The landcover with highest occurrence within the grid
var worldcoverResampled = worldcover
.reduceResolution({
reducer: ee.Reducer.mode(),
maxPixels: 1024
})
.reproject({
crs: gridProjection
});
// Select grids for the following classes
// | Class Name | Value |
// | Forests | 10 |
// | Shrubland | 20 |
// | Grassland | 30 |
// | Cropland | 40 |
// | Mangroves | 95 |
var landCoverMask = worldcoverResampled.eq(10)
.or(worldcoverResampled.eq(20))
.or(worldcoverResampled.eq(30))
.or(worldcoverResampled.eq(40))
.or(worldcoverResampled.eq(95));
var predictedImageMasked = predictedImage
.updateMask(landCoverMask);
Map.addLayer(predictedImageMasked, gediVis, 'Predicted AGBD (Masked)');
หน่วยของค่า GEDI AGBD คือเมกะกรัมต่อเฮกตาร์ (Mg/ha) หากต้องการหา AGB ทั้งหมด เราจะคูณแต่ละพิกเซลด้วยพื้นที่ในหน่วยเฮกตาร์ แล้วนำค่าเหล่านั้นมารวมกัน
var pixelAreaHa = ee.Image.pixelArea().divide(10000);
var predictedAgb = predictedImageMasked.multiply(pixelAreaHa);
var stats = predictedAgb.reduceRegion({
reducer: ee.Reducer.sum(),
geometry: geometry,
scale: gridScale,
maxPixels: 1e10,
tileScale: 16
});
// Result is a dictionary with key for each band
var totalAgb = stats.getNumber('agbd');
print('Total AGB (Mg)', totalAgb);
ลองใช้สคริปต์ฉบับเต็มสำหรับบทแนะนำนี้ในโปรแกรมแก้ไขโค้ดของ Earth Engine