임베딩 필드는 분류에 사용되는 것과 동일한 방식으로 회귀의 특징 입력/예측 변수로 사용할 수 있습니다.
이 튜토리얼에서는 64D 삽입 필드 레이어를 지상 바이오매스 (AGB)를 예측하는 다중 회귀 분석의 입력으로 사용하는 방법을 알아봅니다.
NASA의 Global Ecosystem Dynamics Investigation (GEDI) 미션은 60m 간격으로 30m 공간 해상도에서 지상 트랜섹트를 따라 LIDAR 측정값을 수집합니다. 회귀 모델에서 예측 변수로 사용될 지상 바이오매스 밀도 (AGBD)의 포인트 추정치가 포함된 GEDI L4A Raster Aboveground Biomass Density 데이터 세트를 사용합니다.
지역 선택
관심 영역을 정의하는 것부터 시작해 보겠습니다. 이 튜토리얼에서는 인도 서부 가트 산맥의 한 지역을 선택하고 다각형을 형상 변수로 정의합니다. 또는 코드 편집기의 그리기 도구를 사용하여 관심 영역 주위에 다각형을 그릴 수 있습니다. 이 다각형은 가져오기에서 geometry 변수로 저장됩니다. 또한 식물이 있는 지역을 쉽게 찾을 수 있는 위성 기본 지도도 사용합니다.
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에서 입력 이미지의 모자이크에는 1도 스케일의 WGS84인 기본 투영이 할당됩니다. 이 모자이크는 튜토리얼 후반부에서 집계되고 다시 투영되므로 원래 투영을 유지하는 것이 좋습니다. 타일 중 하나에서 투영 정보를 추출하고 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%)을 모두 삭제합니다.
- Copernicus GLO-30 디지털 고도 모드 (DEM)에 따라 경사가 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 측정값의 수평 정확도는 +/- 9m입니다. 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);
학습 특성 추출
학습 특성을 추출할 입력 데이터가 준비되었습니다. 회귀 모델에서 위성 임베딩 밴드를 종속 변수 (예측 변수)로 사용하고 GEDI AGBD 값을 독립 변수 (예측)로 사용합니다. 각 픽셀에서 일치하는 값을 추출하여 학습 데이터 세트를 준비할 수 있습니다. 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 값과 모델의 예측 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);