স্যাটেলাইট এমবেডিং ডেটাসেটের সাথে রিগ্রেশন

GitHub এ সম্পাদনা করুন
রিপোর্ট সমস্যা
পৃষ্ঠার ইতিহাস
এই টিউটোরিয়ালটি স্যাটেলাইট এমবেডিং ডেটাসেটের টিউটোরিয়ালের একটি সিরিজের অংশ, এছাড়াও ভূমিকা , অ-সুপারভাইজড ক্লাসিফিকেশন , সুপারভাইজড ক্লাসিফিকেশন এবং সিমিলারিটি সার্চ দেখুন।

এম্বেডিং ক্ষেত্রগুলিকে রিগ্রেশনের জন্য বৈশিষ্ট্য ইনপুট/ভবিষ্যদ্বাণী হিসাবে ব্যবহার করা যেতে পারে যেভাবে সেগুলি শ্রেণিবিন্যাসের জন্য ব্যবহৃত হয়।

এই টিউটোরিয়ালে, আমরা শিখব কিভাবে 64D এম্বেডিং ফিল্ড লেয়ারগুলিকে ইনপুট হিসেবে ব্যবহার করতে হয় মাল্টিপল রিগ্রেশন অ্যানালাইসিসের উপর-গ্রাউন্ড বায়োমাস (এজিবি) ভবিষ্যদ্বাণী করে।

NASA-এর গ্লোবাল ইকোসিস্টেম ডায়নামিক্স ইনভেস্টিগেশন (GEDI) মিশন 60 মিটার বিরতিতে 30 মিটার স্থানিক রেজোলিউশনে গ্রাউন্ড ট্রান্সেক্ট বরাবর LIDAR পরিমাপ সংগ্রহ করে। আমরা GEDI L4A রাস্টার অ্যাবোভগ্রাউন্ড বায়োমাস ডেনসিটি ডেটাসেট ব্যবহার করব যেখানে উপরের গ্রাউন্ড বায়োমাস ঘনত্ব (AGBD) এর বিন্দু অনুমান রয়েছে যা রিগ্রেশন মডেলে পূর্বাভাসিত পরিবর্তনশীল হিসাবে ব্যবহার করা হবে।

একটি অঞ্চল নির্বাচন করুন

আসুন আগ্রহের একটি অঞ্চল সংজ্ঞায়িত করে শুরু করি। এই টিউটোরিয়ালের জন্য, আমরা ভারতের পশ্চিমঘাটের একটি অঞ্চল বেছে নেব এবং একটি বহুভুজকে জ্যামিতি পরিবর্তনশীল হিসাবে সংজ্ঞায়িত করব। বিকল্পভাবে, আপনি আগ্রহের অঞ্চলের চারপাশে একটি বহুভুজ আঁকতে কোড এডিটরের অঙ্কন সরঞ্জামগুলি ব্যবহার করতে পারেন যা আমদানিতে জ্যামিতি পরিবর্তনশীল হিসাবে সংরক্ষণ করা হবে। আমরা স্যাটেলাইট বেসম্যাপও ব্যবহার করি, যা গাছপালা অঞ্চলগুলি সনাক্ত করা সহজ করে তোলে।

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 জোনের প্রজেকশনে পরিবেশন করা হয়। ফলস্বরূপ, আমরা একাধিক স্যাটেলাইট এম্বেডিং টাইলস পাই যা আগ্রহের অঞ্চলকে কভার করে। একটি একক চিত্র পেতে, আমাদের তাদের মোজাইক করতে হবে। আর্থ ইঞ্জিনে, ইনপুট চিত্রগুলির একটি মোজাইককে ডিফল্ট প্রজেকশন নির্ধারণ করা হয়, যা 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%)
  • কোপারনিকাস 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 বায়োমাস পর্যবেক্ষণ

নমুনা এবং সমষ্টি ইনপুট

একটি রিগ্রেশন মডেল প্রশিক্ষিত করার জন্য পিক্সেলের নমুনা নেওয়ার আগে, আমরা একই পিক্সেল গ্রিডে ইনপুটগুলিকে পুনরায় নমুনা করি এবং পুনরায় প্রজেক্ট করি। 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 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());

একটি রিগ্রেশন মডেল প্রশিক্ষণ

আমরা এখন মডেল প্রশিক্ষণের জন্য প্রস্তুত. আর্থ ইঞ্জিনের অনেক ক্লাসিফায়ার শ্রেণীবিভাগ এবং রিগ্রেশন উভয় কাজের জন্য ব্যবহার করা যেতে পারে। যেহেতু আমরা একটি সাংখ্যিক মান (একটি শ্রেণীর পরিবর্তে) ভবিষ্যদ্বাণী করতে চাই - আমরা 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)');


চিত্র: মুখোশযুক্ত অ-ভেজিটেড এলাকার সাথে পূর্বাভাসিত AGBD

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

আর্থ ইঞ্জিন কোড এডিটরে এই টিউটোরিয়ালটির জন্য সম্পূর্ণ স্ক্রিপ্টটি চেষ্টা করুন