הודעה: כל הפרויקטים הלא מסחריים שנרשמו לשימוש ב-Earth Engine לפני 15 באפריל 2025 חייבים לעבור אימות של הזכאות לשימוש לא מסחרי כדי לשמור על הגישה ל-Earth Engine.

דף זה תורגם על ידי Cloud Translation API.

רגרסיה עם מערך נתוני הטמעה של לוויין

עריכה ב-GitHub

דיווח על בעיה

היסטוריית הדף

מחברים: spatialthoughts

המדריך הזה הוא חלק מסדרת מדריכים בנושא מערך הנתונים Satellite Embedding. כדאי לעיין גם במדריכים מבוא, סיווג ללא פיקוח, סיווג עם פיקוח וחיפוש דמיון.

אפשר להשתמש בשדות של הטמעה כקלט של תכונות או כמשתנים לחיזוי לרגרסיה, בדיוק כמו שמשתמשים בהם לסיווג.

במדריך הזה נלמד איך להשתמש בשכבות של שדות הטמעה 64D כקלט לניתוח רגרסיה מרובה לחיזוי ביומסה מעל פני הקרקע (AGB).

במסגרת המשימה של נאס"א Global Ecosystem Dynamics Investigation (GEDI), נאספים נתוני LIDAR לאורך חתכי קרקע ברזולוציה מרחבית של 30 מ' במרווחים של 60 מ'. נשתמש במערך הנתונים GEDI L4A Raster Aboveground Biomass Density שמכיל אומדנים נקודתיים של צפיפות ביומסה מעל פני הקרקע (AGBD), שישמשו כמשתנה החזוי במודל הרגרסיה.

בחירת אזור

נתחיל בהגדרת אזור עניין. במדריך הזה נבחר אזור ב-Western Ghats בהודו ונגדיר פוליגון כמשתנה הגיאומטרי. אפשר גם להשתמש בכלי השרטוט בעורך הקוד כדי לשרטט פוליגון סביב האזור שמעניין אתכם. הפוליגון יישמר כמשתנה הגיאומטריה בקטע Imports. אנחנו משתמשים גם במפת בסיס של לוויין, שמאפשרת לאתר בקלות אזורים עם צמחייה.

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 עם קנה מידה של מעלה אחת. בהמשך המדריך נצטרך לצבור ולשנות את ההטלה של הפסיפס הזה, ולכן כדאי לשמור את ההטלה המקורית. אנחנו יכולים לחלץ את פרטי ההטלה מאחת מהמשבצות ולהגדיר אותם בפסיפס באמצעות הפונקציה 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 m. הבעיה הזו מתרחשת כשמנסים להתאים בין ערכי 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));

הקרנה מחדש וצבירה של פיקסלים היא פעולה יקרה, ולכן מומלץ לייצא את התמונה המוערמת שמתקבלת כנכס ולהשתמש בתמונה שחושבה מראש בשלבים הבאים. הפעולה הזו תעזור לכם להתגבר על השגיאות computation timed out או user memory exceeded כשעובדים עם אזורים גדולים.

// 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);

חילוץ תכונות של אימון

הנתונים שלנו מוכנים לחילוץ תכונות האימון. אנחנו משתמשים בפסי Satellite Embedding כמשתנים תלויים (משתני חיזוי) ובערכי 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 מתועדים לעומת ערכים חזויים של המודל

יצירת חיזויים לערכים לא ידועים

אחרי שנגיע למודל שביעות רצון, נוכל להשתמש במודל המאומן כדי ליצור חיזויים במיקומים לא ידועים מתוך התמונה שמכילה פסי חיזוי.

// 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.