رگرسیون با مجموعه داده های جاسازی ماهواره ای

در GitHub ویرایش کنید
گزارش مشکل
تاریخچه صفحه
نویسنده(ها): اندیشه های فضایی
این آموزش بخشی از مجموعه ای از آموزش های مربوط به مجموعه داده های جاسازی ماهواره است، همچنین به مقدمه ، طبقه بندی بدون نظارت ، طبقه بندی نظارت شده و جستجوی مشابه مراجعه کنید.

فیلدهای جاسازی شده را می توان به عنوان ورودی/پیش بینی کننده ویژگی برای رگرسیون به همان روشی که برای طبقه بندی استفاده می شود استفاده کرد.

در این آموزش، ما یاد خواهیم گرفت که چگونه از لایه های میدان جاسازی 64 بعدی به عنوان ورودی برای تحلیل رگرسیون چندگانه برای پیش بینی زیست توده بالای زمین (AGB) استفاده کنیم.

ماموریت جهانی بررسی دینامیک اکوسیستم ناسا (GEDI) اندازه گیری های LIDAR را در امتداد ترانسکت های زمینی با وضوح فضایی 30 متر در فواصل 60 متر جمع آوری می کند. ما از مجموعه داده‌های چگالی زیست توده روی زمین رستر GEDI L4A استفاده خواهیم کرد که حاوی تخمین‌های نقطه‌ای از چگالی زیست توده بالای زمین (AGD) است که به عنوان متغیر پیش‌بینی‌شده در مدل رگرسیون استفاده می‌شود.

یک منطقه را انتخاب کنید

بیایید با تعریف منطقه مورد علاقه شروع کنیم. برای این آموزش، منطقه ای را در گات غربی هند انتخاب می کنیم و چند ضلعی را به عنوان متغیر هندسه تعریف می کنیم. همچنین، می‌توانید از ابزار طراحی در ویرایشگر کد برای ترسیم چند ضلعی در اطراف ناحیه مورد نظر استفاده کنید که به عنوان متغیر هندسه در 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');

مجموعه داده های Satellite Embedding را آماده کنید

تصاویر 64 باندی تعبیه شده ماهواره ای به عنوان پیش بینی کننده رگرسیون استفاده خواهد شد. ما مجموعه داده های Satellite Embedding را بارگذاری می کنیم، تصاویر را برای سال و منطقه انتخاب شده فیلتر می کنیم.

var embeddings = ee.ImageCollection('GOOGLE/SATELLITE_EMBEDDING/V1/ANNUAL');

var embeddingsFiltered = embeddings
  .filter(ee.Filter.date(startDate, endDate))
  .filter(ee.Filter.bounds(geometry));

تصاویر Satellite Embedding در کاشی‌ها شبکه‌بندی می‌شوند و در طرح‌بندی برای مناطق 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 درصد براساس حالت ارتفاعی دیجیتال کوپرنیک GLO-30 (DEM)

در نهایت، تمام اندازه‌گیری‌های باقی‌مانده را برای دوره زمانی و منطقه مورد نظر انتخاب می‌کنیم و یک موزاییک ایجاد می‌کنیم.

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 آماده شده است

نمونه‌برداری مجدد و ورودی‌ها

قبل از نمونه‌برداری از پیکسل‌ها برای آموزش یک مدل رگرسیون، ورودی‌ها را مجدداً نمونه‌برداری می‌کنیم و مجدداً به همان شبکه پیکسلی بازتاب می‌دهیم. اندازه گیری های 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);

استخراج ویژگی های آموزشی

ما داده های ورودی خود را برای استخراج ویژگی های آموزشی آماده کرده ایم. ما از باندهای جاسازی ماهواره ای به عنوان متغیرهای وابسته (پیش بینی کننده) و مقادیر GEDI AGD به عنوان متغیر مستقل (پیش بینی شده) در مدل رگرسیون استفاده می کنیم. ما می توانیم مقادیر همزمان را در هر پیکسل استخراج کنیم و مجموعه داده آموزشی خود را آماده کنیم. تصویر 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);


شکل: مشاهده شده در مقابل مقادیر AGD پیش بینی شده مدل

ایجاد پیش بینی برای مقادیر ناشناخته

هنگامی که از مدل راضی هستیم، می‌توانیم از مدل آموزش‌دیده برای ایجاد پیش‌بینی در مکان‌های ناشناخته از تصویر حاوی باندهای پیش‌بین استفاده کنیم.

// We set the band name of the output image as 'agbd'
var predictedImage = stackedResampled.classify({
  classifier: model,
  outputName: 'agbd'
});

تصویر حاوی مقادیر پیش بینی شده AGD در هر پیکسل اکنون برای صادرات آماده است. در قسمت بعدی از این برای تجسم نتایج استفاده خواهیم کرد.

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


شکل: AGD پیش بینی شده. سبزهای تیره تر نشان دهنده تراکم زیست توده پیش بینی شده بیشتر است

زیست توده کل را تخمین بزنید

ما اکنون مقادیر AGD را برای هر پیکسل تصویر پیش‌بینی کرده‌ایم و می‌توان از آن برای تخمین کل ذخایر زیست توده بالای زمین (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)');


شکل: AGD پیش بینی شده با مناطق غیر پوشش گیاهی پوشانده شده

واحدهای مقادیر GEDI AGD مگاگرم در هکتار (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);

اسکریپت کامل این آموزش را در ویرایشگر کد موتور زمین امتحان کنید .