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

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

התאמה בין Landsat ETM+‎ ל-OLI

עריכה ב-GitHub

דיווח על בעיה

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

Author(s): jdbcode

אזהרה: התהליכים שמתוארים במדריך הזה להרמוניזציה של נתוני Landsat ETM+‎ ו-OLI הם מיושנים ולא מומלצים או נחוצים כשעובדים עם נתוני השתקפות פני השטח של Landsat Collection 2. מידע נוסף זמין בסעיף השאלות הנפוצות בנושא הרמוניזציה של Landsat. הנתונים של אוסף 1 הוצאו משימוש, והקוד במדריך הזה לא יפעל יותר ולא יעודכן.

פתיחה ב-Code Editor

המדריך הזה עוסק בהתאמה של נתוני השתקפות מפני השטח של Landsat ETM+ לנתוני השתקפות מפני השטח של Landsat OLI. הוא מספק:

פונקציית טרנספורמציה ספקטרלית
פונקציות ליצירת נתונים שמוכנים לניתוח
דוגמה להמחשה של סדרת נתונים מבוססת זמן

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

הערה: המקדמים להמרת ETM+‎ ל-OLI חלים גם על TM. לכן, במדריך הזה, כל אזכור של ETM+‎ זהה לאזכור של TM.

מידע על Landsat

‫Landsat היא תוכנית לוויינית להדמיה שאוספת תמונות של כדור הארץ ברזולוציה בינונית מאז 1972. התוכנית הזו היא תוכנית התצפית הארוכה ביותר על כדור הארץ מהחלל, והיא מספקת תיעוד זמני חשוב לזיהוי מגמות מרחביות-זמניות בשינוי הנוף. במדריך הזה נעשה שימוש בנתונים של מכשירים מסוג Thematic Mapper (TM),‏ Enhanced Thematic Mapper Plus (ETM+) ו-Operational Land Imager (OLI). הם קשורים זה לזה באופן הדוק, ויחסית קל לאחד אותם לסדרת זמן עקבית שמפיקה רשומה רציפה משנת 1984 ועד היום, בקצב של 16 ימים לכל חיישן עם רזולוציה מרחבית של 30 מטר. מכשיר הסורק הרב-ספקטרלי (MSS) מרחיב את התיעוד של Landsat אחורה עד 1972, אבל לא נעשה בו שימוש במדריך הזה. הנתונים שלו שונים מאוד, ולכן קשה לשלב אותו עם חיישנים מאוחרים יותר. דוגמאות להרמוניזציה בכל החיישנים אפשר למצוא במאמרים Savage et al. (2018) ו-Vogeler et al. (2018).

למה צריך הרמוניזציה

Roy et al. (2016) demonstrate that there are small, but potentially significant differences between the spectral characteristics of Landsat ETM+ and OLI, depending on application. הסיבות לרצון בהרמוניזציה של מערכי נתונים כוללות: יצירת סדרת זמנים ארוכה שכוללת את Landsat TM,‏ ETM+‎ ו-OLI, יצירת קומפוזיציות תוך-שנתיות קרובות לתאריך כדי לצמצם את ההשפעות של תצפיות חסרות מרווחים של ETM+‎ SLC-off ומיסוך של עננים וצללים, או הגדלת תדירות התצפיות בסדרת זמנים. מידע נוסף זמין בכתב היד המקושר שלמעלה.

הוראות

פונקציות

בהמשך מפורטות סדרת פונקציות שנדרשות כדי להמיר את ETM+‎ ל-OLI ולייצר נתונים שמוכנים לניתוח, שישמשו בקטע 'יישום' במדריך הזה להצגת ההיסטוריה הספקטרלית-זמנית של פיקסל.

הרמוניזציה

ההרמוניזציה מתבצעת באמצעות טרנספורמציה ליניארית של המרחב הספקטרלי של ETM+‎ למרחב הספקטרלי של OLI, בהתאם למקדמים שמוצגים בטבלה 2 של Roy et al. (2016), מקדמי רגרסיית OLS. המקדמים שמתאימים לכל פס מוגדרים במילון הבא עם קבועי התמונה של השיפוע (slopes) ונקודת החיתוך עם הציר (itcps). שימו לב: ערכי נקודת החיתוך עם ציר ה-y מוכפלים ב-10,000 כדי להתאים לקנה המידה של נתוני ההחזרת האור של USGS 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 כוללת פס שמזהה פיקסלים כעננים וצל של עננים, ומגדיר אותם כ-null.

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) מצאו שבקרב 13 אינדקסים/פסים ספקטרליים, ל-NBR היה יחס האות לרעש הגבוה ביותר בהקשר של הפרעות ביערות (אות) ברחבי ארה"ב. הערך מחושב באמצעות הפונקציה הבאה.

function calcNbr(img) {
  return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}

באפליקציה שלכם, אתם יכולים לבחור להשתמש באינדקס ספקטרלי אחר. כאן אפשר לשנות או להוסיף אינדקס נוסף.

שילוב פונקציות

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

// 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), שעבר שינוי כלשהו בשנות ה-80 וה-90, ושריפה בעוצמה גבוהה בשנת 2011.

תחום עניין

איור 1. מיקום ותו באתר לדוגמה של אזור עניין. יער עצי מחט בוגרים בצפון מערב האוקיינוס השקט במדרון הצפוני של הר הוד, אורגון, ארה"ב. תמונות באדיבות: Google Earth,‏ USDA Forest Service,‏ 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 surface reflectance collections עבור 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 כדי לוודא שהפעולה לא חורגת מהמגבלות של מספר הפיקסלים המקסימלי, הזיכרון או הזמן הקצוב לתפוגה. אפשר גם להחליף את פונקציית הצמצום median בסטטיסטיקה המועדפת. מידע נוסף זמין בקטע סטטיסטיקות של אזור בתמונה במדריך למפתחים.

עכשיו אפשר להשתמש בפונקציה ui.Chart.feature.groups כדי לאשר את האוסף. הפונקציה מצפה לקבל אוסף ואת שמות המאפיינים של אובייקט האוסף למיפוי לתרשים. במקרה הזה, המאפיין system:time_start משמש כמשתנה בציר ה-x, והמאפיין NBR משמש כמשתנה בציר ה-y. שימו לב שהמאפיין NBR נקרא על ידי הפונקציה reduceRegion שלמעלה על סמך שם הלהקה בתמונה שעבר צמצום. אם משתמשים בפס אחר (לא 'NBR'), צריך לשנות את הארגומנט של שם המאפיין של ציר ה-y בהתאם. לבסוף, מגדירים את משתנה הקיבוץ (הסדרה) ל-'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, כפי שאפשר לראות ב
עם זאת, חשוב לציין שהשריפה (Dollar Lake) התרחשה בספטמבר 2011. השגיאה בשנה שבה זוהה השינוי נובעת מכך שטווח התאריכים המורכב השנתי הוא מיולי עד אוגוסט. שינויים שמתרחשים אחרי התקופה הזו לא מזוהים עד לתצפית הזמינה הבאה שאינה null, שיכולה להיות בשנה הבאה או מאוחר יותר.
התאוששות הצמחייה (עלייה בתגובת ה-NBR) מתחילה שנתיים אחרי שזוהה אובדן משמעותי של NBR.

פעולות על ציר הזמן – כל התצפיות

איור 2. תרשים של סדרת זמנים של תגובה ספקטרלית לפי פיקסל יחיד עם ייצוג מ-Landsat TM,‏ ETM+‎ ו-OLI. תמונות TM ו-ETM+‎ עוברות הרמוניזציה ל-OLI באמצעות טרנספורמציה לינארית. התמונות הן מיולי עד אוגוסט, והן מסוננות כדי להציג רק תמונות באיכות גבוהה.

תיצור תרשים של סדרות זמנים שמציג את החציון השנתי

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

השלב הראשון הוא לקבץ את התמונות לפי שנה. מוסיפים לכל תמונה מאפיין חדש בשם 'year' על ידי מיפוי הגדרת האוסף 'year' מכל תמונה 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 באמצעות רגרסיה של OLS. פונקציות לכל אפשרויות התרגום זמינות בסקריפט של כלי עריכת הקוד או בסקריפט המקור ב-GitHub.

השריפה באגם דולר

השריפה שמוזכרת במדריך הזה היא השריפה באגם דולר, שהתרחשה באמצע ספטמבר 2011 במדרונות הצפוניים של הר הוד באורגון, ארה"ב. מידע נוסף על התכונה הזו זמין בפוסט הזה בבלוג ובמאמר Varner et al. (2015).