تحذير: إنّ الإجراءات الموضّحة في هذا البرنامج التعليمي لتنسيق بيانات Landsat ETM+ وOLI قديمة ولا يُنصح بها أو لا تكون ضرورية عند العمل مع بيانات انعكاس السطح من Landsat Collection 2. لمزيد من المعلومات، يُرجى الاطّلاع على إدخال الأسئلة الشائعة حول توحيد بيانات Landsat. تم إيقاف بيانات المجموعة 1 نهائيًا، ولن يعمل الرمز البرمجي الوارد في هذا البرنامج التعليمي أو يتم تعديله.
يتناول هذا البرنامج التعليمي عملية مطابقة بيانات انعكاس السطح التي تم جمعها باستخدام جهاز Landsat ETM+ مع بيانات انعكاس السطح التي تم جمعها باستخدام جهاز Landsat OLI. توفّر هذه الخدمة ما يلي:
- دالة تحويل طيفي
- دوال لإنشاء بيانات جاهزة للتحليل
- مثال على تصور السلسلة الزمنية
يهدف هذا الدليل إلى تقديم إرشادات شاملة حول كيفية تنسيق وعرض سلسلة زمنية إقليمية لأكثر من 35 عامًا من بيانات Landsat، ويمكن تطبيقها على الفور على المناطق التي تهمّك.
يُرجى العِلم أنّ معاملات تحويل ETM+ إلى OLI تنطبق على TM أيضًا. وبالتالي، فإنّ الإشارة إلى ETM+ في هذا البرنامج التعليمي هي مرادفة لـ TM.
لمحة عن Landsat
Landsat هو برنامج لتصوير الأقمار الصناعية يجمع صورًا للأرض بدقة متوسطة منذ عام 1972. وباعتباره أطول برنامج لمراقبة الأرض من الفضاء، يقدّم هذا البرنامج سجلاً زمنيًا قيّمًا لتحديد المؤشرات المكانية والزمانية للتغيّرات في المناظر الطبيعية. يتم استخدام بيانات أدوات Thematic Mapper (TM) وEnhanced Thematic Mapper Plus (ETM+) وOperational Land Imager (OLI) في هذا البرنامج التعليمي. وهي مرتبطة ارتباطًا وثيقًا ويسهل نسبيًا دمجها في سلسلة زمنية متسقة تنتج سجلًا مستمرًا من عام 1984 حتى الآن، بمعدل تكرار يبلغ 16 يومًا لكل جهاز استشعار وبدقة مكانية تبلغ 30 مترًا. توسّع أداة "الماسح الضوئي المتعدد الأطياف" نطاق بيانات Landsat ليشمل الفترة منذ 1972، ولكن لا يتم استخدامها في هذا البرنامج التعليمي. وتختلف بياناته تمامًا، ما يجعل عملية الدمج مع أجهزة الاستشعار اللاحقة أمرًا صعبًا. يمكنك الاطّلاع على Savage et al. (2018) وVogeler et al. (2018) للحصول على أمثلة على التوافق بين جميع أجهزة الاستشعار.
أهمية التوافق
أظهرت دراسة أجراها "روي" وآخرون (2016) أنّ هناك اختلافات صغيرة ولكنها قد تكون مهمة بين الخصائص الطيفية لجهاز Landsat ETM+ وجهاز OLI، وذلك حسب التطبيق. تشمل الأسباب التي قد تدفعك إلى تنسيق مجموعات البيانات ما يلي: إنشاء سلسلة زمنية طويلة تغطي بيانات Landsat TM وETM+ وOLI، أو إنشاء مركّبات سنوية داخلية قريبة من التاريخ لتقليل تأثيرات البيانات غير المتوفّرة من فجوات ETM+ SLC وإخفاء السحب والظلال، أو زيادة معدّل تكرار البيانات المرصودة ضمن سلسلة زمنية. يُرجى الاطّلاع على المخطوطة المرتبطة أعلاه للحصول على مزيد من المعلومات.
التعليمات
الدوال
في ما يلي سلسلة من الدوال اللازمة لمواءمة بيانات ETM+ مع بيانات OLI وإنشاء بيانات جاهزة للتحليل سيتم استخدامها في قسم التطبيق من هذا البرنامج التعليمي لعرض السجلّ الطيفي الزمني لبكسل.
التنسيق
يتم تحقيق التوافق من خلال التحويل الخطي للمساحة الطيفية لـ ETM+ إلى المساحة الطيفية لـ OLI وفقًا للمعاملات المقدَّمة في Roy et al. (2016) الجدول 2 لمعاملات الانحدار الخطي العادي.
يتم تحديد معاملات النطاق في القاموس التالي باستخدام ثوابت الصور الخاصة بالميل (slopes) والتقاطع (itcps). يُرجى العِلم أنّه يتم ضرب قيم نقطة التقاطع مع المحور الصادي في 10,000 لتتطابق مع مقياس بيانات انعكاس سطح Landsat التابع لهيئة المسح الجيولوجي الأمريكية.
var coefficients = {
itcps: ee.Image.constant([0.0003, 0.0088, 0.0061, 0.0412, 0.0254, 0.0172])
.multiply(10000),
slopes: ee.Image.constant([0.8474, 0.8483, 0.9047, 0.8462, 0.8937, 0.9071])
};
أسماء نطاقات ETM+ وOLI لنافذة الاستجابة الطيفية نفسها غير متساوية ويجب توحيدها. تختار الدوال التالية نطاقات الانعكاس فقط والنطاق pixel_qa من كل مجموعة بيانات، وتعيد تسميتها وفقًا لنطاق الطول الموجي الذي تمثله.
// Function to get and rename bands of interest from OLI.
function renameOli(img) {
return img.select(
['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']);
}
// Function to get and rename bands of interest from ETM+.
function renameEtm(img) {
return img.select(
['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']);
}
أخيرًا، حدِّد دالة التحويل التي تطبّق النموذج الخطي على بيانات ETM+، وتعيّن نوع البيانات على Int16 لضمان التوافق مع OLI، وتعيد ربط النطاق pixel_qa لاستخدامه لاحقًا في إخفاء السحب والظلال.
function etmToOli(img) {
return img.select(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'])
.multiply(coefficients.slopes)
.add(coefficients.itcps)
.round()
.toShort()
.addBands(img.select('pixel_qa'));
}
إخفاء الغيوم والظلال
يجب إخفاء السُحب وظلالها في البيانات الجاهزة للتحليل. تستخدم الدالة التالية CFmask (Zhu et al., 2015)
pixel_qa، تم تضمين نطاق مع كل صورة انعكاس سطحي من Landsat USGS لضبط
وحدات البكسل التي تم تحديدها على أنّها سحابة وظل سحابة على قيمة فارغة.
function fmask(img) {
var cloudShadowBitMask = 1 << 3;
var cloudsBitMask = 1 << 5;
var qa = img.select('pixel_qa');
var mask = qa.bitwiseAnd(cloudShadowBitMask)
.eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0));
return img.updateMask(mask);
}
احتساب المؤشر الطيفي
يستخدم التطبيق القادم مؤشر الطيف "نسبة الاحتراق العادية" (NBR) لتمثيل السجلّ الطيفي لبكسل في غابة تضرّرت بسبب حريق غابات. يتم استخدام نسبة NBR لأنّ Cohen et al. (2018) وجدوا أنّ نسبة NBR، من بين 13 مؤشرًا/نطاقًا طيفيًا، حققت أعلى نسبة إشارة إلى الضوضاء فيما يتعلق باضطراب الغابات (الإشارة) في جميع أنحاء الولايات المتحدة. ويتم احتسابها من خلال الدالة التالية.
function calcNbr(img) {
return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
في تطبيقك الخاص، يمكنك اختيار استخدام مؤشر طيفي مختلف. يمكنك هنا تعديل الفهرس أو إضافة فهرس آخر.
دمج الدوال
تحديد دالة تضمين لكل مجموعة بيانات تجمع كل الدوال المذكورة أعلاه لتسهيل تطبيقها على مجموعات الصور الخاصة بها
// Define function to prepare OLI images.
function prepOli(img) {
var orig = img;
img = renameOli(img);
img = fmask(img);
img = calcNbr(img);
return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}
// Define function to prepare ETM+ images.
function prepEtm(img) {
var orig = img;
img = renameEtm(img);
img = fmask(img);
img = etmToOli(img);
img = calcNbr(img);
return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}
في تطبيقك، يمكنك اختيار تضمين وظائف أو استبعادها. عدِّل هذه الدوال حسب الحاجة.
مثال على السلسلة الزمنية
يمكن ربط دالتَي التغليف prepOli وprepEtm أعلاه بمجموعات انعكاس السطح من Landsat لإنشاء بيانات جاهزة للتحليل من أجهزة استشعار متعددة بهدف عرض التسلسل الزمني الطيفي لوحدة بكسل أو منطقة من وحدات البكسل. في هذا المثال، ستنشئ سلسلة زمنية تزيد عن 35 عامًا وتعرض السجلّ الطيفي لبكسل واحد. تتعلّق هذه البكسل بتاريخ حديث لمساحة من غابات الصنوبريات الناضجة في شمال غرب المحيط الهادئ (الشكل 1) التي تعرّضت لبعض الاضطرابات في الثمانينيات والتسعينيات وحريق كبير في عام 2011.
الشكل 1. الموقع الجغرافي وطابع الموقع لمنطقة الاهتمام النموذجية غابة صنوبرية ناضجة في شمال غرب المحيط الهادئ على المنحدر الشمالي لجبل هود، أوريغون، الولايات المتحدة الأمريكية الصور مقدَّمة من: Google Earth وخدمة الغابات التابعة لوزارة الزراعة الأمريكية وLandsat وCopernicus.
تحديد منطقة اهتمام (AOI)
ونتيجة هذا التطبيق هي سلسلة زمنية من بيانات Landsat المرصودة لبكسل معيّن. يتم استخدام عنصر ee.Geometry.Point لتحديد موقع البكسل.
var aoi = ee.Geometry.Point([-121.70938, 45.43185]);
في تطبيقك، يمكنك اختيار بكسل مختلف. يمكنك إجراء ذلك من خلال استبدال إحداثيات خطوط الطول والعرض المذكورة أعلاه. بدلاً من ذلك، يمكنك تلخيص السجلّ الطيفي لمجموعة من وحدات البكسل باستخدام تعريفات أخرى للكائنات ee.Geometry، مثل ee.Geometry.Polygon() وee.Geometry.Rectangle().
يُرجى الاطّلاع على قسم
الأشكال الهندسية
في "دليل المطوّر" للحصول على مزيد من المعلومات.
استرداد مجموعات بيانات أجهزة استشعار Landsat
يمكنك الحصول على مجموعات انعكاس السطح من Landsat USGS لأجهزة الاستشعار OLI وETM+ وTM.
انتقِل إلى الروابط لمعرفة المزيد عن كل مجموعة بيانات.
var oliCol = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
var etmCol = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');
var tmCol = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR');
تحديد فلتر لمجموعة صور
حدِّد فلترًا لتقييد مجموعات الصور حسب الحدود المكانية للمنطقة محل الاهتمام وموسم ذروة التمثيل الضوئي والجودة.
var colFilter = ee.Filter.and(
ee.Filter.bounds(aoi), ee.Filter.calendarRange(182, 244, 'day_of_year'),
ee.Filter.lt('CLOUD_COVER', 50), ee.Filter.lt('GEOMETRIC_RMSE_MODEL', 10),
ee.Filter.or(
ee.Filter.eq('IMAGE_QUALITY', 9),
ee.Filter.eq('IMAGE_QUALITY_OLI', 9)));
يمكنك تغيير معايير الفلترة هذه كما تريد. حدِّد خصائص أخرى للفلترة من خلال علامة التبويب "خصائص الصور" في صفحات وصف البيانات المرتبطة أعلاه.
إعداد المجموعات
ادمج المجموعات وطبِّق الفلتر واربط الدالة prepImg بجميع الصور. نتيجة المقتطف التالي هي صورة ee.ImageCollection واحدة
تحتوي على صور من مجموعات أجهزة الاستشعار OLI وETM+ وTM التي تستوفي
معايير الفلترة، ويتم معالجتها لتصبح جاهزة للتحليل باستخدام مؤشر NBR.
// Filter collections and prepare them for merging.
oliCol = oliCol.filter(colFilter).map(prepOli);
etmCol = etmCol.filter(colFilter).map(prepEtm);
tmCol = tmCol.filter(colFilter).map(prepEtm);
// Merge the collections.
var col = oliCol.merge(etmCol).merge(tmCol);
إنشاء رسم بياني لسلسلة زمنية يعرض جميع الملاحظات
يمكن تقييم التوافق بين أجهزة الاستشعار بشكل نوعي من خلال رسم جميع الملاحظات كمخطط بياني مبعثر يتم فيه تمييز جهاز الاستشعار حسب اللون. توفّر دالة
ui.Chart.feature.groups في Earth Engine هذه الأداة، ولكن يجب أولاً إضافة قيمة NBR المرصودة إلى كل صورة كسمة. اربط دالة تقليل المنطقة بمجموعة الصور التي تحسب متوسط جميع وحدات البكسل المتقاطعة مع منطقة الاهتمام وتضبط النتيجة كخاصية للصورة.
var allObs = col.map(function(img) {
var obs = img.reduceRegion(
{geometry: aoi, reducer: ee.Reducer.median(), scale: 30});
return img.set('NBR', obs.get('NBR'));
});
في تطبيقك، إذا كانت منطقة الاهتمام كبيرة، يمكنك زيادة قيمة scale وتحديد وسيطات bestEffort وmaxPixels وtileScale للدالة reduceRegion لضمان عدم تجاوز العملية للحدود القصوى لوحدات البكسل أو الذاكرة أو المهلة. يمكنك أيضًا استبدال أداة تقليل الوسيط بالإحصاءات المفضّلة لديك. راجِع قسم
إحصاءات منطقة صورة
في دليل المطوِّر لمزيد من المعلومات.
يمكن الآن قبول المجموعة من خلال الدالة ui.Chart.feature.groups.
تتوقّع الدالة مجموعة وأسماء خصائص عنصر المجموعة ليتم ربطها بالرسم البياني. في هذا المثال، تعمل السمة "system:time_start" كمتغير للمحور x، بينما تعمل السمة "NBR" كمتغير للمحور y. يُرجى العِلم أنّ الدالة reduceRegion أعلاه قد سمّت السمة NBR استنادًا إلى اسم النطاق في الصورة الذي تم تقليله. إذا كنت تستخدم نطاقًا مختلفًا (ليس "NBR")، غيِّر وسيطة اسم السمة الخاصة بالمحور الصادي وفقًا لذلك. وأخيرًا، اضبط متغيّر التجميع (السلسلة) على "SATELLITE"، والذي يتم تعيين ألوان مميزة له.
var chartAllObs =
ui.Chart.feature.groups(allObs, 'system:time_start', 'NBR', 'SATELLITE')
.setChartType('ScatterChart')
.setSeriesNames(['TM', 'ETM+', 'OLI'])
.setOptions({
title: 'All Observations',
colors: ['f8766d', '00ba38', '619cff'],
hAxis: {title: 'Date'},
vAxis: {title: 'NBR'},
pointSize: 6,
dataOpacity: 0.5
});
print(chartAllObs);
سيظهر رسم بياني مشابه للشكل 2 في وحدة التحكّم بعد بعض الوقت من المعالجة. يُرجى ملاحظة ما يلي:
- تتوفّر عدة ملاحظات في السنة.
- تتألف السلسلة الزمنية من ثلاثة أجهزة استشعار مختلفة.
- لا توجد انحيازات واضحة في أجهزة الاستشعار.
- يختلف معدّل تكرار الملاحظات سنويًا.
- هناك بعض التباين خلال العام، وهو متسق على مدار السنوات ومحدود نسبيًا. عند تجربة ذلك باستخدام مؤشر NDVI (img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');) ، يظهر تفاوت أكبر بكثير في الاستجابة داخل السنة وبين السنوات.
- يظلّ معدّل NBR مرتفعًا لمعظم السلاسل الزمنية مع حدوث بعض الاضطرابات الطفيفة.
- يؤدي حريق الغابات إلى انخفاض كبير وسريع في استجابة NBR، كما يتضح في
- يُرجى العِلم أنّ الحريق (بحيرة دولار) وقع في سبتمبر 2011. يعود الخطأ في سنة الرصد إلى أنّ النطاق الزمني السنوي المركّب يمتد من تموز (يوليو) إلى آب (أغسطس)، ولا يتم رصد التغييرات التي تحدث بعد هذه الفترة إلا عند توفّر الملاحظة التالية غير الفارغة، والتي قد تكون في العام التالي أو بعده.
- يبدأ تعافي الغطاء النباتي (زيادة استجابة مؤشر نسبة العجز الطبيعي) بعد عامين من رصد فقدان كبير في مؤشر نسبة العجز الطبيعي.
الشكل 2. مخطط السلسلة الزمنية للاستجابة الطيفية لوحدة بكسل واحدة مع تمثيل من Landsat TM وETM+ وOLI يتم تنسيق صور TM وETM+ مع OLI من خلال تحويل خطي. الصور من تموز (يوليو) إلى آب (أغسطس) وتمت فلترتها للحصول على صور عالية الجودة.
إنشاء رسم بياني لسلسلة زمنية يعرض المتوسط السنوي
لتبسيط السلسلة الزمنية وإزالة التشويش منها، يمكن تطبيق عملية تقليل الملاحظات خلال السنة. في هذه الحالة، يتم استخدام الوسيط.
الخطوة الأولى هي تجميع الصور حسب السنة. أضِف خاصية "السنة" جديدة إلى كل صورة من خلال الربط بإعداد المجموعة "السنة" من ee.Image.Date لكل صورة.
var col = col.map(function(img) {
return img.set('year', img.date().get('year'));
});
يمكن استخدام السمة الجديدة "السنة" لتجميع الصور من السنة نفسها من خلال عملية ربط. يؤدي الربط بين مجموعة صور مميزة من سنة معيّنة والمجموعة الكاملة إلى توفير قوائم من السنة نفسها، ويمكن إجراء عمليات تقليل المتوسطات منها. قسِّم المجموعة الكاملة إلى مجموعة من ممثلي السنوات المختلفة.
var distinctYearCol = col.distinct('year');
حدِّد فلترًا وعملية ربط لتحديد السنوات المتطابقة بين مجموعة السنوات المميزة (distinctYearCol) والمجموعة الكاملة (col). سيتم حفظ التطابقات في سمة باسم "year_matches".
var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');
طبِّق عملية الربط وحوِّل FeatureCollection الناتج إلى ImageCollection.
var joinCol = ee.ImageCollection(join.apply(distinctYearCol, col, filter));
// Apply median reduction among matching year collections.
var medianComp = joinCol.map(function(img) {
var yearCol = ee.ImageCollection.fromImages(img.get('year_matches'));
return yearCol.reduce(ee.Reducer.median())
.set('system:time_start', img.date().update({month: 8, day: 1}));
});
غيِّر ee.Reducer أعلاه إذا كنت تفضّل استخدام إحصاء آخر غير الوسيط لتمثيل مجموعة الصور في سنة معيّنة.
أخيرًا، أنشئ رسمًا بيانيًا يعرض قيم NBR المتوسطة من مجموعات الملاحظات الصيفية السنوية. غيِّر إحصاءات تقليل المنطقة حسب ما يناسبك.
var chartMedianComp = ui.Chart.image
.series({
imageCollection: medianComp,
region: aoi,
reducer: ee.Reducer.median(),
scale: 30,
xProperty: 'system:time_start',
})
.setSeriesNames(['NBR Median'])
.setOptions({
title: 'Intra-annual Median',
colors: ['619cff'],
hAxis: {title: 'Date'},
vAxis: {title: 'NBR'},
lineWidth: 6
});
print(chartMedianComp);
سيظهر رسم بياني مشابه للشكل 3 في وحدة التحكّم بعد بعض الوقت من المعالجة.
الشكل 3. مخطط السلسلة الزمنية لمتوسط الاستجابة الطيفية في الصيف لوحدة بكسل واحدة مع تمثيل من Landsat TM وETM+ وOLI يتم تنسيق صور TM وETM+ مع OLI من خلال تحويل خطي. الصور من تموز (يوليو) إلى آب (أغسطس) وتمت فلترتها للحصول على صور عالية الجودة.
معلومات إضافية
دوال التحويل البديلة
يقدّم الجدول 2 في Roy et al. (2016) معاملات الانحدار بطريقة المربعات الصغرى العادية (OLS) وطريقة الانحدار بتقليل المدى (RMA) لتحويل انعكاس سطح ETM+ إلى انعكاس سطح OLI والعكس. يوضّح البرنامج التعليمي أعلاه عملية تحويل ETM+ إلى OLI فقط باستخدام الانحدار الخطي البسيط. يمكنك العثور على دوال جميع خيارات الترجمة في نص "محرّر الرموز" البرمجية أو نص المصدر على GitHub.
حريق بحيرة الدولار
الحريق المذكور في هذا البرنامج التعليمي هو حريق Dollar Lake الذي اندلع في منتصف أيلول (سبتمبر) 2011 على المنحدرات الشمالية لجبل Hood في ولاية أوريغون، الولايات المتحدة الأمريكية. يمكنك الاطّلاع على منشور المدوّنة هذا وVarner et al. (2015) للحصول على مزيد من المعلومات حول هذا الموضوع.