吐血整理如何在Google Earth Engine上寫迴圈 五個程式碼例項詳細拆解

绕shoe三匝發表於2024-07-18

在這裡同步一篇本人的原創文章。原文釋出於2023年釋出在知乎專欄,轉移過來時略有修改。全文共計3萬餘字,希望幫助到GEE小白快速進階。

引言

這篇文章主要解答GEE中.map().iterate()函式的用法。

首先解答一個疑問,為什麼需要自己寫迴圈?確實,GEE 為各種資料型別提供了無數常用的內建函式,對這些方法做排列組合足以應對大多數使用場景,演算法效率也頗佳。在有內建函式的情況下,我永遠推薦內建函式。

然而,對策趕不上需求,總有時候需要部署個性化的演算法。舉個例子,資料聚合問題:現有一個小時級的氣溫資料集 hourlyTemp_imgcol,每小時整點提供一張近地表氣溫影像,要求借此生成一個每日平均均氣溫資料集 dailyTem_imgcol,每天提供一張日均近地表氣溫影像。那麼,GEE 有沒有這樣的內建函式:

var dailyTemp_imgcol = hourlyTemp_imgcol.aggregate('daily');

讓我們一行程式碼出結果呢?沒有。

以上的例子是細分學科下的任務,聽起來還是挺稀罕的。那我就再舉一個程式設計初學者一定會接觸的演算法,數列遞推:由一個初始值按照一定的遞推公式,生成一個固定長度的列表。GEE有沒有類似這樣的遞推函式:

var List = Number.RECURSIVE_SEQUENCE(function, list_length);

來一步到位實現這個簡單任務呢?很遺憾,也是沒有的。

因此,儘管 GEE 內建函式庫豐富,學會自己部署迴圈仍然很有必要。這是跨過初學者階段、獨立開展專案的必經之路。

在第一個例子中,我們需要按時間迴圈,每個日期計算一幅平均氣溫影像。在第二個例子中,我們需要按列表索引位置迴圈,每個索引位置上計算一個遞推值。

第一個例子代表了一類問題,處理函式作用在每個元素上,讓每個元素被對映到一個新的東西上,這樣的話一定是多個輸入和多個輸出,並且有一一對應的關係。關於這種迴圈,我們使用.map()操作來處理,請閱讀本文第(一)節。

第二個例子代表了另一類問題,處理函式的任務是產生累積的效果,後面步驟的執行會依賴前面步驟的產出。關於這一種迴圈,我們使用.iterate()操作來處理,請閱讀本文第(二)節。

有些初學者會說,哪裡需要那麼麻煩,我以不變應萬變,老師教的for迴圈就是無腦易上手,再說GEE又不是不支援。我想說,這種習慣在技術上確實可行,但是如果是放在 GEE 上部署,會嚴重拖慢計算效率。在展示如何寫 GEE 式迴圈之前,我專門開一個第(〇)節,探討一下這個問題。非常建議初學者閱讀,會幫助你更好地理解,GEE 如何處理使用者提交的任務。

然而並不是說任何情況都不能使用for迴圈。 有些情況下,因為.map().iterate()的技術限制,迴圈必須在本地端進行。我會舉一下這種應用場景的例子,放在第(三)節。

讓我們開始吧。

(〇) 為什麼GEE不推薦for、while迴圈?

GEE 常被描述為“地理雲端計算平臺”。那麼,之所以被稱為雲端計算,就是因為操作和計算是在伺服器上發生的,而不是在自己的電腦(本地)上完成的。向伺服器告知計算任務時所用的“語言”,稱為API。目前已有JavaScript API(即Code Editor)、Python API、Julia API。

-- 那麼老哥,forwhile是跟GEE伺服器對話的語言嗎?

-- 不是。

-- 難道說GEE伺服器不看for語句?

-- 它不看。

-- 唉不對啊,我明明在Code Editor上面寫了for迴圈卻還能執行的,你說伺服器不看for,難不成還能是瀏覽器給我跑的?

-- 恭喜你答對了,還真是瀏覽器給跑的。

要解釋這個問題,需要對GEE的工作機制有一個大概的瞭解。

作為小白,我們在Code Editor上寫出幾行程式碼並點選Run的一刻,可能很少想過這之後發生了什麼。實際上,你的程式碼其實並不是直接抄一份送給伺服器去讀,而是要先在本地端瀏覽器上去做解析,重寫成一套伺服器聽得懂的請求,然後才把請求送上去。

本地翻譯的過程,在於解析每一個由使用者宣告的變數,把它們各自整理成一套能指導伺服器算出結果的結構——你可以理解成是每一個雲端物件在本地端的代理,是影子。它們會跟雲端的、承載實際資料和計算的物件相對應。

在《Deferred Execution(延遲執行)》的頁面上[1],谷歌解釋了這一本地解析的過程:

When you write a script in Earth Engine (either JavaScript or Python), that code does NOT run directly on Earth Engine servers at Google. Instead, the client library encodes the script into a set of JSON objects, sends the objects to Google and waits for a response. Each object represents a set of operations required to get a particular output, an image to display in the client, for example.
當您在Earth Engine上編寫指令碼時(無論JavaScript還是Python),該程式碼不會直接在谷歌的Earth Engine伺服器上執行。相反,客戶端庫將指令碼編碼為一組JSON物件,將物件傳送給谷歌並等待響應。每個物件表示獲得特定輸出所需的一組操作,例如,要在客戶端上顯示的影像。

進一步解釋一下:程式碼寫好之後,只會在本地做解析。本地端解析之後,會生成一些JSON物件。這些JSON物件是做什麼用的呢?是客戶端寫給伺服器的字條:“如此如此算一遍,把結果發給我。”

由此可以看出,本地端解析程式碼所能做的,只是解析程式碼,幫伺服器提前理解該怎麼做事而已,不涉及任何實際資料。它編寫的JSON檔案,只是雲端資料在本地端的影子。生成這個影子的目的,或者說這個影子裡所容納的資訊,就是指揮Google伺服器以某某運算組合依次操作那個雲端上的物件。

這樣就解釋了為什麼for迴圈只在本地端有意義。我們在for迴圈裡寫的一切,是一種過程而不是一種物件,只會平添本地端的解析負擔。for迴圈對伺服器上計算過程產生影響之前,是要本地客戶端來徒手拆迴圈的。

說句題外話,這其實很像在編纂紀傳體史書。歷史的原始記錄就像我們寫的程式碼一樣,關注的是事情從頭到尾怎麼進行,是一種過程。而紀傳體史書則關注主要人物個體如何發展,這更像GEE伺服器的關注點,是一種物件。把流水賬般的事項記錄(過程),編製成人物個體的發展歷程(物件及其操作順序),並對互相影響的人物各自立傳、一起成書,這些就是史官(和GEE API)的工作。

理解了本地端和伺服器端物件的區別和聯絡,讓我們再來討論以for為代表的本地端迴圈,為什麼不被推薦。

舉個例子,我們這裡有一個重複1萬次的for迴圈,裡面的內容是讓一個ee.Number(伺服器端的數)每次加 i:

// 語法沒錯但跑不通的例子,只是因為用了 for 迴圈
var num = ee.Number(0);
for (var i=1; i<=10000; i++) {
  num = num.add(i);
}
print(num);

已知伺服器不看也看不懂for迴圈,那麼這個拆迴圈的工作,一定是在本地完成的。瀏覽器需要認真的把1萬次迭代走完。然而,它在迭代時可不是替伺服器執行實際計算,而是在字條裡寫一步請求:嘿伺服器,這一步,你執行加 i 操作。於是它把這句話又寫了9999次。在迴圈結束時,它出的活將是一個超長的、一萬個處理步驟的字條。

這就是瀏覽器上的客戶端API在面對for迴圈時負擔的壓力。瀏覽器表示寫字條太苦逼了,求求你還是讓我來算數吧。

看著是不是替自己的瀏覽器捏把汗呢?實際上,這個程式碼在我的電腦上根本跑不通,因為這個預計超長的JSON物件,是超出了我的瀏覽器寫字條的能力的(報錯:Maximum call stack size exceeded)。

一旦選擇去寫for迴圈,那麼制約程式碼執行效率的因素,就不在於伺服器的計算速度了,而是本地端的翻譯能力。而這還不是最差的情況。最不推薦的做法是for迴圈中還帶有.getInfo()函式,向伺服器請求發回這一步計算結果,它讓本該只是影子的本地變數取回它在伺服器端對應的值。這個函式會打斷非同步計算的工作方式,在當前值被髮送回本地前,一切其它步驟的計算都不能接著進行,同時伺服器也閒置了。如果被套進很多層for迴圈中,使用者實際的感受就是,本地等待時間非常長,瀏覽器無響應(甚至提示崩潰)。

反面案例:

// for 迴圈與 getInfo 合用則是更大的錯誤!
// 以下程式碼能跑通,但親測可以讓瀏覽器卡死數十分鐘
// 不要執行!!!
var num = ee.Number(0);
for (var i=1; i<=10000; i++) {
  num = ee.Number(num).add(i).getInfo();
}
print(num);

那麼,使用.map().iterate()就能繞過這個問題嗎?是的。這是因為,它們兩個都是伺服器端的函式,是伺服器可以直接讀懂的東西。這樣一來,呼叫它們也就無需讓客戶端去拆迴圈,只需一次解析,剩下的事情就交給伺服器了。瀏覽器壓力大減。

要是forwhile迴圈哪怕尚有一成的用處,我的這篇文章也就沒有什麼必要了。但我可以拍著胸脯說,除了極少數情況(見第三節),都是可以尋找到對本地迴圈的替代方案的。這裡推薦閱讀GEE的官方教程文件《Functional Programming Concepts[2](函數語言程式設計的概念),教你把習慣的程序式程式設計做法(for迴圈、if/else條件語句、累積迭代)重構為適合GEE的程式碼。

(一) 逐項對映用.map()

引例

有多個型別的 ee 物件都帶有.map()操作:ee.List(列表),ee.Dictionary(字典),ee.ImageCollection(影像集),ee.FeatureCollection(要素集)。借用Python語言中的概念,它們都是可迭代物件iterables。ee.List和ee.Dictionary類似於Python的List和Dictionary。ee.ImageCollection和ee.FeatureCollection則非常類似於Python的集合,set。

巧的是,在Python中也有一個map函式處理可迭代物件,用法如下,是一個讓數列逐項加1的指令碼。

""" map(fun, iter) 
 \* fun - 處理函式,map將每個元素傳遞給它
 \* iter - 可迭代物件,map作用於它的每個元素上
"""
list1 = [1, 2, 3, 4]
list2 = list( map(lambda x: x+1, list1) )
list2     # [2, 3, 4, 5]

這個指令碼,放在 GEE Code Editor 中的話,應該這樣寫:

/** ee.List.map(baseAlgorithm)
 * this: list - 輸入物件,此例中為列表,map作用於它的每個元素上
 * baseAlgorithm - 處理函式,map將每個元素傳遞給它 
*/
var list1 = ee.List([1, 2, 3, 4]);
var list2 = list1.map(
  function(x){ 
    return ee.Number(x).add(1);
  }
);
print(list2); // [2, 3, 4, 5]

可以看到,map.map()在兩個平臺上有著一樣的本質功能,都是逐個元素處理,逐一對映。也就是說,把處理函式作用在輸入物件的每個元素上,並且在輸出物件中建立一一對應的新元素。

二者在編寫結構上有一些細節上的差別:

  1. 儘管都是需要兩個輸入,Python的map是在它的括號裡套起來處理函式和輸入物件。而GEE中的.map() 則只套處理函式,跟在輸入物件後。
  2. Python的map可以同時傳遞進來多個可迭代物件(以上示例中沒有展示),而GEE的.map()則只能處理一個物件。
  3. Python使用lambda語句建立匿名函式,而GEE Code Editor用function語句。由於採用function語句,GEE Code Editor使用者可以在.map()的括號裡寫一個帶更多步驟的匿名函式(以上示例中沒有展示)。
  4. Python的map的直接輸出是一個map object,需要用container套起來去賦予型別,比如寫成 list(map(...)) 。而GEE中的.map()則返回一個與輸入物件型別相同的物件,進來是什麼型別出去就還是什麼型別。

這種處理是並行化進行的。在上面的例子中,我們並不關心是數字1先被處理,還是數字4先,只在乎它們被處理完後有沒有被擺在正確的位置。4個數字或許是同時被各自送入不同的 CPU 核,同時被輸出,同時被填放。又或許是先處理3個,再處理1個,齊活之後再一起上菜。又或許是先處理2個再處理2個……具體怎麼排程的,那是後臺的事,對我們使用者其實毫無影響。

但值得放在腦子中的資訊是,這種處理不會是像for迴圈那樣,一個一個依次操作,而是以並行化的方式提高效率。記住這一點,後面要考的。

簡單的.map()操作常用於對影像集中的影像、要素集中的向量要素作批次處理,編輯、處理後,生成對應的新資料。在更復雜的情況中,通常是自定義的特殊處理任務,我們需要以給定的時間、空間、或關鍵字等資訊為依據,對其它資料作聚合,如統計日平均氣溫,即氣溫資料在時間維度的聚合。對於這兩種.map()操作的應用場景,這裡給出例1和例2兩個假想的資料處理任務以供參考。

例1 GIS資料的批次編輯

這裡給出的例子是對於影像集(ee.ImageCollection)進行批次編輯。對要素集(ee.FeatureCollection,或者稱向量集)的批次編輯雖然在任務性質上會有差異,但在對.map()的使用方式上同理,所以就不予贅述了。

我們設想這樣一種任務,對影像集裡的每一個影像,都進行一次掩膜操作,把不想要的像元篩除掉。這樣的任務常見於定量遙感,常常是資料準備階段中的一步預處理。掩膜圖層常常是特定的土地覆蓋型別(如森林、水體、農田),或氣候條件閾值(如溫度、降水量),或有云/無雲的位置。我這裡舉兩個非常最常用的Landsat預處理操作:1. 使用質量控制波段,對Landsat 8地表反射率資料集進行掩膜。2. 使用官方給定的放縮乘數(scale factor)和偏置(offset),將資料集裡的源資料轉換成能夠使用的地表反射率。

程式碼連結

首先,取得Landsat 8地表反射率資料集(注意是Collection 2)。這裡僅取同一個地點上2個月的影像,一共4張,因為取多了處理時間會拉長。4張影像的雲量分別為:89.8%,0.05%,21.39%,75.44%。

/////////////////////////////////////////////////////////////
// 獲取Landsat 8 Collection 2 Tier 1地表反射率(Level 2)影像集
// 獲取到的影像ID及雲量:
// LANDSAT/LC08/C02/T1_L2/LC08_012031_20210711 - 89.8%
// LANDSAT/LC08/C02/T1_L2/LC08_012031_20210727 - 0.05%
// LANDSAT/LC08/C02/T1_L2/LC08_012031_20210812 - 21.39%
// LANDSAT/LC08/C02/T1_L2/LC08_012031_20210828 - 75.44%
var landsat8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
                .filter(ee.Filter.and(
                  ee.Filter.eq('WRS_PATH', 12),
                  ee.Filter.eq('WRS_ROW', 31)))
                .filterDate('2021-07-01', '2021-09-01');

其次,寫一個函式來實現:對一幅輸入進來的Landsat 8影像,利用它自身的質量控制波段,去除雲和雲影。質量控制波段是一幅數值影像。在二進位制化後,它的每一個(或2個)數位指示像元質量的某一個屬性。我們需要第1至第4數位,分別是dilated cloud(檢測到的雲範圍的放大範圍)、cirrus cloud(捲雲)、cloud(明顯的雲)、cloud shadow(雲影)。

/////////////////////////////////////////////////////////////
////////////// 雲掩模函式 //////////////
function prep_l8 (l8_image) {
  // 轉換成影像,以保證資料型別正確
  l8_image = ee.Image(l8_image);
  // 使用QA數位1-4,構建掩模圖層(取1是不受雲影響的像元,取0是受雲影響的像元)
  // Bit 0: fill (0)
  // Bit 1: dilated cloud 
  // Bit 2: cirrus cloud
  // Bit 3: cloud
  // Bit 4: cloud shadow
  var l8_qa = l8_image.select("QA_PIXEL"); 
  var dilated_cloud_mask = l8_qa.bitwiseAnd(1 << 1).eq(0);
  var cirrus_cloud_mask = l8_qa.bitwiseAnd(1 << 2).eq(0);
  var cloud_mask = l8_qa.bitwiseAnd(1 << 3).eq(0);
  var cloud_shadow_mask = l8_qa.bitwiseAnd(1 << 4).eq(0);
  var mask = dilated_cloud_mask
              .and(cirrus_cloud_mask)
              .and(cloud_mask)
              .and(cloud_shadow_mask);
  // 放縮和偏置,並恢復在乘法加法中丟失的後設資料(metadata)
  var scale = 0.0000275;
  var offset = -0.2;
  var l8_image_scaled = ee.Image(l8_image.multiply(scale).add(offset)
                                  .copyProperties(l8_image))
                            .set(
                              "system:id", l8_image.get("system:id"),
                              "system:version", l8_image.get("system:version")
                              );
  // 掩模並返回影像
  var l8_image_scaled_masked = l8_image_scaled.updateMask(mask);
  return l8_image_scaled_masked;
}

接下來執行批次編輯,也就是讓上面的處理函式作用在影像集的每個影像身上。具體的寫法,是讓.map()括起來處理函式的名字(或者可以在.map()裡寫函式),再跟隨在待處理的影像集後面即可。

/////////////////////////////////////////////////////////////
// 批次預處理,顯示結果
landsat8 = landsat8.map(prep_l8);

var visualization = {
  bands: ['SR_B4', 'SR_B3', 'SR_B2'],
  min: 0.0,
  max: 0.1,
};
Map.addLayer( ee.Image(landsat8.toList(10).get(0)), visualization, "2021-07-11" );
Map.addLayer( ee.Image(landsat8.toList(10).get(1)), visualization, "2021-07-27" );
Map.addLayer( ee.Image(landsat8.toList(10).get(2)), visualization, "2021-08-12" );
Map.addLayer( ee.Image(landsat8.toList(10).get(3)), visualization, "2021-08-28" );

顯示這4張經過預處理後的影像。可以發現,我們用一步操作為4張影像去除掉了受雲影響的像元,可以說.map()的老本行是GIS資料批次處理實在不為過。

例2 時間維度上的資料聚合

GEE 上能實現的聚合不僅限於求和。做乘積、求平均、求方差、找極值等等,只要用得上就都能做。我們這裡以求平均為例子說明做資料聚合的技巧。

在時間維度上做資料聚合,輸入到演算法中資訊有兩個:時間集、資料集。不同於在批次編輯操作中的做法,我們這回不能再讓.map()當資料集的“跟屁蟲”,而是要轉而去跟隨時間集。

原因很好理解。.map()操作,輸入輸出的維度必須是相同的。我們如果想實現一個時間出一張圖,那就必須讓時間成為.map()的輸入。(至於資料集,它可以去當處理函式的輸入。這裡有個非常好用的技巧我接下來會在處理函式程式碼處細說)

做法是把時間維度放到一個 ee.List 列表物件中。因為 ee.List 也帶有.map()操作,我們的程式設計靈活性有了保證。

考慮到ee.List.map()的輸出變數是另一個 ee.List,我們有必要對結果進行變換——重建所需的資料型別,比如轉換成 ee.ImageCollection 或 ee.FeatureCollection。

這裡展示的程式碼將實現我們在引言中提到的,將逐小時氣溫資料轉換為每日平均氣溫。可以點選此連結在 GEE Code Editor 中檢視。

首先,定義輸入量:

  1. 逐小時近地表氣溫資料,由實時中尺度分析(the Real-Time Mesoscale Analysis,RTMA)資料集[3]來提供。這是由美國國家氣象局(National Weather Service,NWS)提供的氣象再分析產品,覆蓋美國大陸48州,空間解析度2500米,時間解析度1小時。
  2. 時間範圍,2022年1月1日至1月3日。初始化一個 ee.List 存放所有日期的起始時刻,以美國東部時間為時區。
var hourly_temp = ee.ImageCollection("NOAA/NWS/RTMA").select("TMP");
var dates = ee.List([
  ee.Date('2022-01-01T00:00:00', 'EST5EDT'),
  ee.Date('2022-01-02T00:00:00', 'EST5EDT'),
  ee.Date('2022-01-03T00:00:00', 'EST5EDT')
  ]);

然後,定義處理函式。在每一個日期上,這個函式要實現三個步驟:首先對 hourly_temp 影像集做篩選,然後對篩選後的影像集作聚合生成一張均值影像,最後重新命名均值影像。

/** 
 * 本函式構造一個內函式並返回
 * 使內函式滿足ee.List.map()對處理函式格式的要求 
 */
function calc_daily_mean(imgcol) { 
  /** 在內函式中完成處理 */
  return function(idate) {
    idate = ee.Date(idate);
    /** 第一步:對影像集按日期做篩選 */
    var filtered_imgcol = imgcol.filterDate(idate);
    /** 第二步:對篩選後的影像集作聚合,生成一張均值影像 */
    var mean_img = filtered_imgcol.mean();
    /** 第三步01:提取原影像集的波段名稱 */
    var band_names = imgcol.first().bandNames();
    /** 第三步02:用“原波段名稱+年月日”重新命名均值影像的波段,如:"TMP20220101" \*/
    band_names = band_names.map(function(name){return ee.String(name).cat(idate.format('YYYYMMdd','EST5EDT'))});
    return mean_img.rename(band_names);
  };
}

終於可以把前面賣的關子揭開了。可以看到,處理函式的結構是,一個外函式返回一個內函式。為什麼要這樣設計呢?我們在引例中提到,GEE 的.map()則只能處理一個輸入變數。這樣就產生了矛盾,因為我們的資料聚合任務必須含有兩個變數:時間集和資料集。當我們已經把時間集指定為.map()的輸入時,資料集似乎就無處可去了。然而,假如我們稍微改造一下處理函式,寫成這樣的外函式套內函式、並返回內函式的結構,就可以做到讓資料集輸入到外函式,時間集裡的日期輸入到內函式。這樣一來,可以同時滿足兩方面看似矛盾的要求。

使用這一處理函式的程式碼如下。第一輸入——時間集 dates ——放在.map()的左側。第二輸入——資料集hourly_temp ——就讓外函式括起來。然後,在整個語句的最外層,用ee.ImageCollection()套起來,把.map()的輸出結果(原本是ee.List)轉換成影像集。

var daily_temp = ee.ImageCollection(dates.map(calc_daily_mean(hourly_temp)));

最後,讓我們在地圖檢視中檢視一下算好的平均氣溫。另外,再生成一個顏色條方便閱讀資料。

var vis = {
  min: -30,
  max: 20,
  palette: ['001137','01abab','e7eb05','620500']
};
Map.addLayer(daily_temp.first(), vis, "2022-01-01 mean temperature");
Map.addLayer(ee.Image(daily_temp.toList(3).get(1)), vis, "2022-01-02 mean temperature");
Map.addLayer(ee.Image(daily_temp.toList(3).get(2)), vis, "2022-01-03 mean temperature");

var legendTitle = ui.Label({
  value: 'Mean temperature (C)',
  style: {fontWeight: 'bold'}
});
var legendLabels = ui.Panel({
  widgets: [
    ui.Label(vis.min, {margin: '4px 8px'}),
    ui.Label(
        ((vis.max-vis.min) / 2+vis.min),
        {margin: '4px 8px', textAlign: 'center', stretch: 'horizontal'}),
    ui.Label(vis.max, {margin: '4px 8px'})
  ],
  layout: ui.Panel.Layout.flow('horizontal')
});
var colorBar = ui.Thumbnail({
  image: ee.Image.pixelLonLat().select(0),
  params: {
    min: 0,
    max: 1,
    palette: vis.palette,
    bbox: [0, 0, 1, 0.1],
  },
  style: {stretch: 'horizontal', margin: '0px 8px', maxHeight: '24px'},
});
Map.add(ui.Panel([legendTitle, colorBar, legendLabels]));


美國大陸部分2022年1月1日平均氣溫地圖。


2022年1月2日


2022年1月3日

.map()易錯點

在GEE上使用.map()時,有必要牢記一些易錯點。.map()裡執行出的錯,報錯反饋可能不會指出問題出在了程式碼第幾行,也許也無法提供足夠精確的診斷資訊(這可能是 GEE 最難用的地方,不過 GEE 似乎正在改進這一點)。把易錯點爛熟於心,至少對於debug是很有幫助的。

首先,輸入物件和元素類別的問題。.map()的輸入變數型別只能是 ee.List、ee.Dictionary、ee.ImageCollection、ee.FeatureCollection 四者之一。除了這些ee物件以外,不能對其它型別使用.map()。並且,輸入物件中的元素經由.map()送入處理函式中時,一律都會被伺服器理解成一個“沒有型別的”雲端物件 ee.ComputedObject,而不是原本的型別,所以在處理函式中需要為元素重新定義物件類別

比如,引例中 ee.List 的元素,在使用.add()進行處理之前,需要被定義成 ee.Number。即使我們把輸入列表中的物件都寫成 ee.Number,在被送入處理函式後也仍然會“失去”型別——見下方的錯誤程式及其報錯。

/** 這是一個錯誤示例 */
/** 說明數值元素在送入處理函式時沒有被理解為ee.Number,而是需要重新定義 */
var list1 = ee.List([ee.Number(1), ee.Number(2), ee.Number(3), ee.Number(4)]);
var list2 = list1.map(
  function(x){ 
    return x.add(1); // 正確寫法:return ee.Number(x).add(1)
  }
);
print(list2); // 報錯:"x.add is not a function"

其次,本地端操作的問題。由於.map()迴圈是完全交給伺服器執行的,它的處理函式里面不能寫printiffor這些程式設計師喜聞樂見的語句。一旦寫上了就會彈出報錯"A mapped function's arguments cannot be used in client-side operations"(對映函式的引數不能在客戶端操作中使用)。print的問題實在解決不了,但是iffor還是有替代品的。例如ee.Algorithms.If可以實現if,不過更推薦的做法是用.filter()操作提前篩選好資料。要是手癢癢想寫for的話,那就疊一個.map().iterate()迴圈好咯。

.map()小結

  1. .map()適用於逐個元素處理、逐一對映的任務。
  2. .map()在 GEE 後臺採用並行化的執行方式。
  3. .map()的兩種常用場景:GIS資料的批次編輯,以及在時間、空間、或關鍵字等維度上進行資訊聚合。
  4. 使用時應注意變數的資料型別,以及避免在處理函式中使用本地端操作。

(二) 累積迭代用.iterate()

前一小節中,我們討論了.map()操作的適用範圍——對映。細心的同學們可能已經發現了,.map()對元素的處理是一種各自獨立的訪問方式,順序無所謂誰先誰後,只要處理好後放置順序正確即可。因此,它的處理函式沒有“步”的概念,也就不能產生按步累積的效果。

但是我們學習for迴圈時,還有一個很重要應用就是按步求累積啊,這可怎麼辦?

這就是.iterate()的功能了。.iterate()操作適用於那些需要考慮狀態逐步積累、逐步改變的計算任務。

這裡提供一個簡單的判斷標準,決定你該不該用.iterate()。如果你的問題必須拆解成子步驟逐步計算,且後步必須依賴前步的計算結果,那麼它就是.iterate()的應對範圍。反之,則優先嚐試使用.map()

適合用.map()去解決的問題較多,而非得用.iterate()的問題較少。.iterate()執行起來是完全按照規定好的1234步來做的。能用.map()去並行執行的任務非要寫成.iterate(),倒也不是出不了結果,但屬實是有點一核有難多核圍觀了。

帶有.iterate()操作的物件型別有3個:ee.List,ee.ImageCollection,ee.FeatureCollection,比.map()少一個 ee.Dictionary。

呼叫.iterate()的程式碼格式比.map()繁瑣一些。 .iterate()需要兩個輸入:處理函式 + 第一步迴圈開始前的初始狀態,二者都是必須的輸入量,後者是.map()沒有的。同樣繁瑣的還有.iterate()對處理函式的要求:迴圈控制變數 + 前一輪的迭代狀態(作為本輪迭代的基準)。

下文中給出用.iterate()實現做累積迭代的一些任務。例3是對儲存在 ee.FeatureCollection 中的要素屬性進行求和的例子。例4匯流演算法是一個相當複雜的迴圈迭代過程,套在大一些的研究區上很容易超出GEE計算限制,因此不建議在GEE上部署。

例3 簡單的迭代求和

這裡舉一個官方文件[4]中的例子。構造一個要素集(ee.FeatureCollection),其中每個要素以感興趣區ROI為幾何形狀,以一個月份的月降水量為欄位,共12個要素。.iterate()承擔的任務是以迭代的方式對12個月降水量求和,得到感興趣區內的年降水量。

https://developers.google.com/earth-engine/apidocs/ee-featurecollection-iterate

首先,資料來源和感興趣區。資料來源選擇了一個按月更新的氣候資料集(影像集ee.ImageCollection),取一年期。感興趣區選在美國新墨西哥州的一片區域。

 /**
 * CAUTION: ee.FeatureCollection.iterate can be less efficient than alternative
 * solutions implemented using ee.FeatureCollection.map or by converting feature
 * properties to an ee.Array object and using ee.Array.slice and
 * ee.Array.arrayAccum methods. Avoid ee.FeatureCollection.iterate if possible.
 */

// Monthly precipitation accumulation for 2020.
var climate = ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE')
                  .filterDate('2020-01-01', '2021-01-01')
                  .select('pr');

// Region of interest: north central New Mexico, USA.
var roi = ee.Geometry.BBox(-107.19, 35.27, -104.56, 36.83);

然後,使用分割槽統計的方式逐一構造向量要素來儲存月降水量,合成一個要素集。這裡使用了.map()操作來實現對每一張月降水量影像出產當月的向量要素。

// A FeatureCollection of mean monthly precipitation accumulation for the
// region of interest.
var meanPrecipTs = climate.map(function(image) {
  var meanPrecip = image.reduceRegion(
      {reducer: ee.Reducer.mean(), geometry: roi, scale: 5000});
  return ee.Feature(roi, meanPrecip)
      .set('system:time_start', image.get('system:time_start'));
});

以上都是準備階段,下面進入正題。首先建立一個處理函式,名為cumsum,以當前要素(變數名currentFeature)為迴圈控制變數,以已經參與過月降水量加和的要素所組成的列表(變數名featureList,型別為ee.List)為前一輪迭代狀態。處理函式所做的事情,是把已有的月降水量的和(featureList最尾端要素的“pr_cumsum”欄位)和當前月份的降水量(currentFeature的“pr”欄位)分別取出並求和,再新建一個要素賦予這個求和值,並把新建的要素新增到featureList列表中。

// A cumulative sum function to apply to each feature in the
// precipitation FeatureCollection. The first input is the current feature and
// the second is a list of features that accumulates at each step of the
// iteration. The function fetches the last feature in the feature list, gets
// the cumulative precipitation sum value from it, and adds it to the current
// feature's precipitation value. The new cumulative precipitation sum is set
// as a property of the current feature, which is appended to the feature list
// that is passed onto the next step of the iteration.
var cumsum = function(currentFeature, featureList) {
  featureList = ee.List(featureList);
  var previousSum = ee.Feature(featureList.get(-1)).getNumber('pr_cumsum');
  var currentVal = ee.Feature(currentFeature).getNumber('pr');
  var currentSum = previousSum.add(currentVal);
  return featureList.add(currentFeature.set('pr_cumsum', currentSum));
};

為了使用這一處理函式,我們以上文中構造好的月降水量要素集為輸入,並且再構造一個迴圈起始時的要素列表(變數名first,將“pr_cumsum”欄位設定為0)。在迴圈結束後,卸磨殺驢,去除掉這個佔位用的要素。

// Use "iterate" to cumulatively sum monthly precipitation over the year with
// the above defined "cumsum" function. Note that the feature list used in the
// "cumsum" function is initialized as the "first" variable. It includes a
// temporary feature with the "pr_cumsum" property set to 0; this feature is
// filtered out of the final FeatureCollection.
var first = ee.List([ee.Feature(null, {pr_cumsum: 0, first: true})]);
var precipCumSum =
    ee.FeatureCollection(ee.List(meanPrecipTs.iterate(cumsum, first)))
        .filter(ee.Filter.notNull(['pr']));

由此我們可以畫出這樣的年內降水量曲線。藍線表示月降水量,紅線表示累積降水量。這樣一來,我們藉助.iterate(),給for迴圈的迭代求和功能找到了替代。

在以上例子中我們利用.iterate()以一種迭代的方式去對所有要素按一個欄位求和,然而這並不是唯一的實現方式。實際上,如果只想要對欄位求和(但不要求同時給出逐時累積曲線)的話,GEE有內建函式.aggregate_sum(property)提供相同的功能,不必要自己編寫迭代程式碼。

實際上,很多我們直覺上認為需要寫累積迴圈的情況,GEE都貼心的提供了內建函式。內建函式不光呼叫更方便,而且往往執行更快。

例4 匯流問題

注:GEE上已入庫全球90米級別的MERIT Hydro資料集[5]。如果不追求更高解析度,建議直接使用現有產品,無需親自執行匯流演算法。

匯流問題地形分析中的經典老問題,但也是一個需要巨量迭代、狂吃記憶體和時間的老大難技術性問題。

匯流問題通常是這樣的:已知一片感興趣區的DEM影像,要生成一張新影像,使得每個像元的值取該位置所對應的集水區(catchment area,或catchment basin)的面積。集水區的定義:對於地表上任何一個地點而言,凡是落在它鄰近某個區域內的雨水,經過不斷匯聚和流動都會流到這個地點,這個雨水降落和匯流的區域就稱為該地點的“集水區”。落在這個集水區邊界以外的雨水,不論怎樣流動,都不會經過這個地點[6]


集水區的概念圖(圖源:Nazhad et al., 2017 [7])。任意一個地點(圖上小紅點和大紅點)都可以有各自對應的集水區(黃圈和綠圈)。一般來說,上游地點擁有的集水區面積較小,下游地點擁有的集水區面積更大。分析集水區有助於推斷河流的形成。

那麼這個匯流問題,有什麼難點呢?

思考一下簡化版的問題,為一個像元,我們稱為A,圈出它的集水區,不能圈多不能圈少。這等同於需要標記全場其它像元BCDEFG…,哪些能給像元A送水,哪些不能。

似乎,是不是要遍歷全場像元呀?是的,而且不止一遍。首先,對任意像元BCDEFG…,要規劃水落在該位置上可能的流向。為了簡化,我們假設水只會選擇向8-鄰接像元流動,並且只會選擇其中的1個像元為流動方向,也就是下坡坡度最大方向的那個像元。所以在第一次遍歷中,我們需要在全圖範圍上,逐像元按照鄰域地形確定下坡坡向。這張坡向圖告訴我們的是,水落下來後緊接著的下一步,只有一步哦,的流動方向。

有了坡向圖/流向圖,能順水推舟圈集水區了嗎?可以了。思路是,在A的鄰域裡逆向查詢流向指向A的像元,標記下來。再對新標記的所有像元,執行同樣的查詢過程,再標記新一批像元,以此類推。直到標出範圍外,或者標無可標(延伸到了流域邊界),則停止。停止這個逆流追溯過程後,把已標記的全部像元轉換成像元總數(或總面積),賦值給A。

所以,為解決一個像元相應集水區的簡化版問題,我們要開動一次追溯。在這個過程中,全場像元被遍歷了幾遍,取決於客觀上有多少個像元給A送水。只要上一步的追溯還能發現新像元,哪怕一個,卷積核就得給我再接著滑。


匯流問題的示意圖,DEM -> 流向圖 -> 集水區。此例子中為所有像元都統計集水區。圖源: [8]

上文所述,是為一個像元圈出集水區的思路,它只是為了方便向大家說明演算法需要考慮的基本問題。事實上,為全場像元都統計集水區,並不需要逐個像元跑追溯,那樣會產生太多重複步驟,浪費計算資源。最佳化的方式,是把下游向上遊“追溯”的思路,切換成上游向下遊“傳播”的思路。第一次迭代中,每個像元向處在自己流向/坡向的唯一鄰居像元,告知自己的存在(或自己的面積),而接收了訊號的像元需要對訊號求和。第二步及之後的每次迭代中,還是向同樣的下游鄰居像元通訊,但告知內容變成,自己在上一步(且僅在上一步)中收集到的、來自更上游像元的通訊總數(或總面積)。這樣,來自上游的“訊號”逐漸向下遊傳遞且被累加起來,每個像元都能在一個迭代步中更新當前已經發現的上游像元總數(或總面積),且不會出現重複計數。直到所有最上游像元的訊號都傳遞到各自相應的最下游,結束迴圈。

但是,技術上減少計算量的手段,並不能改變匯流問題需要反覆執行迭代操作的本質。科學家針對這一問題提出過很多種演算法以求最佳化這一問題的解決方案,但該問題上仍然缺乏簡單高效的對策。在對範圍廣大、解析度高的DEM資料做處理時,匯流演算法的處理時間能以十小時甚至百小時為單位。

這裡姑且把“傳播”演算法的實現方式給貼出來。純粹是為了技術展示,並不是真的號召大家在GEE上寫這種東西。

匯流演算法的GEE程式碼來自GIS Stack Exchange論壇上的問答[9][10],作者為使用者Bert Coerver,轉載程式碼請註明原始出處。

首先,作者手動新增了一個感興趣區,位於印度。

var geometry1 = ee.Geometry.Polygon(
        [[[81.32080078125, 19.114029215761995],
          [80.958251953125, 18.999802829053262],
          [80.947265625, 18.672267305093758],
          [81.71630859375, 18.396230138028827],
          [82.09808349609375, 18.747708976598744],
          [82.02667236328125, 19.072501451715087]]]);

這裡我需要註明一下,匯流演算法的作用區域,一般要至少覆蓋一整個流域,否則集水區面積的計算結果是會偏小的。這裡程式碼作者隨手畫的這個範圍應該是沒有覆蓋流域,因此計算結果僅供教學,可能不具備現實中的指示作用。


感興趣區示意圖,由程式碼作者Bert Coerver標畫

其次,定義輸入量:

  1. 迭代次數預設為120,這個值要不小於區域內最長的那條匯水路徑的長度(像元數),以保證所有傳播都能完成。
  2. 集水區的計量方式採用了像元數(單位:個),讀者可以自行採用程式碼註釋內提供的替代方式,把計量方式改成面積、單位改成平方千米。
  3. 初始化了一張名為data的影像,用作上文所述的自上而下傳播演算法的初始迭代量。
  4. 坡向圖採用GEE Data Catalog中現成的HydroSHEDS水流流向資料集(30弧秒解析度,約900米)[11],像元值標示水流在該位置的流向:1=東,2=東南,4=南,8=西南,16=西,32=西北,64=北,128=東北。
// Set the number of routing iterations:
var iteration_steps = ee.Number(120);

// Give an image with some data to rout:
var data = ee.Image(1);
// or E.G.: data = ee.Image.pixelArea().divide(1e6);

// Give the bandname from 'data' to rout:
var bandname = ee.String("constant");
// or E.G.: var bandname = ee.String("area");

// Give unit of the data:
var unit = "pixels";
// or E.G.: unit = "km2";

// Give image with flow directions:
var dir = ee.Image("WWF/HydroSHEDS/30DIR").select(["b1"]);

// Set Area-Of-Interest:
var geom = geometry1;

// Set output data type:
var typ = "uint32";

接下來,重新命名迭代量data影像的波段為“rout”(可能是單詞route/routed的簡寫),該波段的初值為1,即1個像元,使用者可以自行在上方程式碼中改用像元面積為計量。複製rout波段,建立一個新波段賦給data影像,重新命名為“summed”,這樣summed波段初值也為1(或自身像元面積)。

rout波段可以理解成每個像元在當前步驟持有的“訊號”,它或來自自身初始化或來自上游傳入,並且這個訊號將在後續步驟中傳遞出去。summed波段可以理解為,每個像元已發現的匯流像元數量或面積(自身也算,如果不想算自身的話,summed可以賦初值為0)。

// Rename the band.
data = data.select([bandname],["rout"]);

// Add a band to store the summation of the routs.
data = data.addBands(data.select(["rout"],["summed"]));

// Set property indicating the current iteration step.
data = data.set("iter_idx", 0);

借用ee.FeatureCollection建立表格(向量資料的屬性表)的能力,建立一個集合來定義方向。這裡面每個向量要素都不含地理參考資訊(geometry欄位為null),因此是“偽要素”,我們專注於屬性表裡的資訊即可。

// Create a feature-collection that describes to what directions (e.g. 128 is NE).
var col = ee.FeatureCollection(ee.List([
          ee.Feature(null, {"weight": ee.List([[0,0,0],[1,0,0],[0,0,0]]),
                            "direction": "E",
                            "dir_id": 1,}),
          ee.Feature(null, {"weight": ee.List([[1,0,0],[0,0,0],[0,0,0]]),
                            "direction": "SE",
                            "dir_id": 2,}),
          ee.Feature(null, {"weight": ee.List([[0,1,0],[0,0,0],[0,0,0]]),
                            "direction": "S",
                            "dir_id": 4,}),
          ee.Feature(null, {"weight": ee.List([[0,0,1],[0,0,0],[0,0,0]]),
                            "direction": "SW",
                            "dir_id": 8,}),
          ee.Feature(null, {"weight": ee.List([[0,0,0],[0,0,1],[0,0,0]]),
                            "direction": "W",
                            "dir_id": 16,}),
          ee.Feature(null, {"weight": ee.List([[0,0,0],[0,0,0],[0,0,1]]),
                            "direction": "NW",
                            "dir_id": 32,}),
          ee.Feature(null, {"weight": ee.List([[0,0,0],[0,0,0],[0,1,0]]),
                            "direction": "N",
                            "dir_id": 64,}),
          ee.Feature(null, {"weight": ee.List([[0,0,0],[0,0,0],[1,0,0]]),
                            "direction": "NE",
                            "dir_id": 128,})]));

⬆可以看到,“direction”和“dir_id”欄位,正是對照著HydroSHEDS水流流向資料集對於流向/坡向的定義,但“weight”欄位中的矩陣則是反過來的。比如,[[0,0,0], [1,0,0], [0,0,0]]應該指向西,卻對應於“東”。這樣設計是因為,上游像元向下遊像元傳播訊號的過程,實際上是靠下游像元的主動收集來實現的(影像卷積操作)。上游像元向東傳送的訊號,需要下游像元從西面去接收。“direction”和“dir_id”欄位所定義的是訊號傳送的去向,而“weight”欄位定義的則是訊號接收的來向。有必要說明一下這個問題以免產生疑惑。

最核心的部分來了:傳播演算法的執行函式。每一步迭代時,正是這個函式來把“訊號”向下遊移動一步。

// The iteration algorithm.
var iterate_route = function(iter_step, data){
    
    // The function to do one routing iteration for one direction.
    var route = function(ft){
  
        // Create the kernel for the current flow direction.
        var kernel = ee.Kernel.fixed({width: 3, 
                                      height: 3, 
                                      weights: ft.get("weight")
                                     });
    
        // Get the number corresponding to the flow direction map.
        var dir_id = ee.Number(ft.get("dir_id"));
        
        // Mask irrelevent pixels.
        var routed = ee.Image(data).select("rout").updateMask(dir.eq(dir_id));
    
        // Move all the pixels one step in the direction currently under consideration.
        var routed_dir = routed.reduceNeighborhood({reducer: ee.Reducer.sum(), 
                                                    kernel: kernel,
                                                    skipMasked: false
                                                   });
                                                     
        return routed_dir;
        };
        
    // Loop over the eight flow directions and sum all the routed pixels together.
    var step = ee.ImageCollection(col.map(route)).sum().rename("rout")
                                                       .reproject(dir.projection())
                                                       .set("iter_idx", iter_step)
                                                       .clip(geom);
    
    // Sum the newest routed pixels with previous ones.
    var summed = step.select("rout").add(ee.Image(data).select("summed")).rename("summed");
    
    // Add the 'rout' and 'summed' bands together in one image.
    var data_next_step = step.addBands(summed);
    
    return data_next_step;
    };

⬆可以看到,這個大函式iterate_route裡巢狀了一個小函式route。小函式的輸入是八方向之一,作用是實現所有像元集體朝一個方向收聽一次,判斷該方向上是否是一個上游像元,收集上游像元帶多少訊號。讓小函式對我們的八方向集合中的每個方向逐一處理(用.map()實現),並對結果求和,也就實現了訊號的一次從上游到下游的傳播。由傳播的結果重建rout波段。傳播完成之後,把rout波段和summed波段求和,用來更新summed波段。

推演一下,對於 i 取任意值,第 i 步迭代之後,來自訊號都向下游流動一步,彙集到了下游像元的rout和summed波段中。而沒有訊號流入的像元(可能是因為自身處在最上游,或是因為自身的上游鄰居在第 i-1 步沒有更上游的訊號流入),則在第 i 步之後rout值為0(這意味著第 i+1 步迭代時無法向下遊傳遞訊號),summed值不變。

定義好執行函式後,就只需呼叫.iterate()執行即可。建立一個名為steps的數列,取值1至120。以steps規定的步驟,按步執行執行函式iterate_route。將輸出影像的資料型別轉換成32位整型。

// Create dictionary to cast output to the chosen datatype.
var cast = ee.Dictionary({"rout" : typ,
                          "summed": typ});

// Create the list with iteration step numbers.
var steps = ee.List.sequence(1, iteration_steps);

// Do the actual iterations.
var full = ee.Image(steps.iterate(iterate_route, data)).cast(cast);

最後,在地圖視窗上展示匯流演算法的輸出結果。

// Add Layer to Map
Map.addLayer(full.select("summed"), {min: 0, max: 1000}, "summed");
Map.centerObject(full.select("summed"))


匯流演算法給出的各像元集水區,計量方式:像元個數。黑色像元約為1-100,灰色像元值數百,亮白色像元1000以上。

後續:在GEE上執行匯流演算法,尤其需要注意預設一個合適的迭代次數。因為GEE的.iterate()函式不存在類似breakcontinue那樣跳過或提前中斷迴圈的語句,所以迭代次數就成為了關乎任務規模的唯一控制量。如果預設的太少,傳播沒有完成就結束了迴圈,集水區計算就會偏小。如果預設的太大,傳播完成後的那些迭代次數就成了沒有必要的空轉,嚴重浪費時間。建議提前估算出全流域內最長的一條匯水路徑的長度,並將迭代次數設定為略大於這條路徑上的像元數量。並且,在演算法的輸出裡,建議加上標記傳播當前狀態的rout波段,作為演算法結束後校驗傳播是否已經客觀完成的參考依據。

另外,需要提醒注意,執行時間過長、記憶體佔用過大的任務可能無法在GEE上順利執行。執行匯流演算法時,流域範圍不宜過大、空間解析度不宜過高。雖然GEE沒有公佈過給每個使用者的記憶體配額,或是執行時間的上限(經驗值是12小時),但GEE官方開設的開發者討論區經常有彙報相關報錯資訊的帖子。對於那些需要執行大流域、高解析度匯流演算法的朋友們來說,GEE也許不是非常穩妥的選擇。如果90米解析度可以接受的話,GEE已經入庫的MERIT Hydro資料集[5]可以直接拿來使用,該資料提供像元數量和麵積兩種集水區計量方式,並且覆蓋全球。

.iterate()小結:

  1. .iterate()適用於需要考慮狀態逐步積累、逐步改變的問題。在不需要按步執行的迴圈問題上,應該優先使用.map()去實現。
  2. 在很多常見的應用場景下,GEE上都存在對.iterate()的替代函式。優先使用GEE內建函式以提高計算效率。
  3. .iterate()的輸入量有兩個:執行函式 + 第一步迴圈開始前的初始狀態。
  4. .iterate()要求執行函式的輸入量為:迴圈控制變數 + 前一輪迭代狀態。輸出量為本輪迭代狀態。
  5. 小心選擇任務規模。提交過於龐大的迭代任務可能會導致超出GEE的記憶體配額或執行時間限制。可以說GEE的這些限制對計算量龐大的迴圈迭代很不友好。如有這方面需求,請謹慎考慮使用GEE。
  6. .iterate()不支援跳過或中斷迴圈,預估任務所需的迭代次數並設計輸入變數的元素數量很重要。
  7. .map()的易錯點,比如在處理函式中需要為元素重新定義物件類別、不支援本地端操作,對.iterate()同樣適用。

(三) 使用本地迴圈的情況

在第(〇)節中,我曾苦口婆心地勸大家不要用for迴圈,甚至為此不惜冒著被大家左滑退出的風險談了些枯燥的、技術性的東西。然而,我必須承認,有的時候本地迴圈就是不得不用。

我們在第(一)節說過,.map()不支援printiffor這些語句,這些限制對於.iteration()也存在。這裡再補充兩條,它們倆內部同樣不支援Export,也不支援Map.addLayer()。前者是將計算結果輸出到Google Drive、GEE Assets、以及Google Cloud Storage,後者是在地圖皮膚上新增資料圖層。我個人覺得,在迴圈裡用Map.addLayer()以便新增圖層,以及用Export以便批次輸出,還是比較常見的需求。既然GEE的迴圈函式不支援,我們就不得不轉而尋求本地迴圈。

在例5中,用一個類似for迴圈的語句forEach來實現本地迴圈,使每個迴圈執行一套獨立的計算,將結果表格輸出到Google Drive中。

例5 批次輸出

本例來自於GIS Stack Exchange上的問答[12],轉載程式碼請註明原始出處。本文對部分原始程式碼進行了一點小修改以避免檔名格式報錯。

先介紹一下這段程式碼的大意。有2個空間位置(83.1499°E, 23.14985°N 和 79.477034°E, 11.015846°N),以及2個時段(1980-2006 和 2024-2060)。任務是對每個空間位置、每個時段,從某個氣候資料庫(NASA NEX-GDDP)中整理出來每日時間序列,各生成1個表格輸出。共計2x2=4個表格。

2個空間位置的座標存放在sites列表中,2個時段的始末年份存放在yearRanges列表中,這些列表都是本地列表(而非ee.List)。

構造了一個自定義函式exportTimeseries來按座標和時段提取資料記錄,函式輸入為一組座標和時段,函式內呼叫Export.table.toDrive來輸出表格。

使用forEach搭配匿名函式,編寫了一個二重迴圈,以便對sites和yearRanges實現迴圈遍歷。在二重迴圈內部執行自定義的exportTimeseries函式。

var sites = [
  [83.1499, 23.14985],
  [79.477034, 11.015846],
];

var yearRanges = [
  [1980, 2006],
  [2040, 2060]
];

sites.forEach(function (site) {
  yearRanges.forEach(function (yearRange) {
    exportTimeseries(site, yearRange);
  });
});

function exportTimeseries(site, yearRange) {
  var startDate = ee.Date.fromYMD(yearRange[0], 1, 1);
  var endDate = ee.Date.fromYMD(yearRange[1], 1, 1);
  
  // get the dataset between date range and extract band on interest
  var dataset = ee.ImageCollection('NASA/NEX-GDDP')
  .filter(ee.Filter.date(startDate,endDate));
  var maximumAirTemperature = dataset.select('tasmax');
  
  // get projection information
  var proj = maximumAirTemperature.first().projection();
  
  var point = ee.Geometry.Point(site);
  
  // calculate number of days to map and extract data for
  var n = endDate.difference(startDate,'day').subtract(1);
  
  // map over each date and extract all climate model values
  var timeseries = ee.FeatureCollection(
    ee.List.sequence(0,n).map(function(i){
      var t1 = startDate.advance(i,'day');
      var t2 = t1.advance(1,'day');
      var feature = ee.Feature(point);
      var dailyColl = maximumAirTemperature.filterDate(t1, t2);
      var dailyImg = dailyColl.toBands();
      // rename bands to handle different names by date
      var bands = dailyImg.bandNames();
      var renamed = bands.map(function(b){
        var split = ee.String(b).split('_');
        return ee.String(split.get(0)).cat('_').cat(ee.String(split.get(1)));
      });
      // extract the data for the day and add time information
      var dict = dailyImg.rename(renamed).reduceRegion({
        reducer: ee.Reducer.mean(),
        geometry: point,
        scale: proj.nominalScale()
      }).combine(
        ee.Dictionary({'system:time_start':t1.millis(),'isodate':t1.format('YYYY-MM-dd')})
      );
      return ee.Feature(point,dict);
    })
  );

  var name = yearRange.join('-') + '_' + [site[0].toPrecision(2),site[1].toPrecision(2)].join('-');
  Map.addLayer(point, null, name);
  Map.centerObject(point,6);
  
  
  // export feature collection to CSV
  Export.table.toDrive({
    collection: timeseries,
    description: 'a_hist_tmax_' + name,
    fileFormat: 'CSV',
  });
}

(四)總結

我用5個例項展示了一些GEE上常見的迴圈的用法。.map()2個,.iterate()2個,本地迴圈1個。它們當然不可能覆蓋全部的場景,不過我想也足夠提供一定的技術性參考了。寫這些主要是為了方便新手查閱,畢竟初學GEE時實在是有著相當陡峭的學習曲線。次要一點的動力就是,希望我的技術分享能給這個逐漸走向內容凋敝的中文網際網路環境一點新鮮血液吧。

寫迴圈可能是衡量GEE使用者是否是熟練工的標準。祝GEE初學者們儘快走過新手期,掌握這個必要技巧。不過,我有必要在結尾重複一下引言中的那句話,在有內建函式的情況下,我永遠推薦內建函式

歡迎與我交流探討!

參考

[1] How Earth Engine Works: Deferred Execution | GEE Guides https://developers.google.com/earth-engine/guides/deferred_execution

[2] Functional Programming Concepts - Google Earth Engine https://developers.google.com/earth-engine/tutorials/tutorial_js_03

[3] RTMA: Real-Time Mesoscale Analysis - Earth Engine Data Catalog https://developers.google.com/earth-engine/datasets/catalog/NOAA_NWS_RTMA

[4] Google Earth Engine - ee.FeatureCollection.iterate https://developers.google.com/earth-engine/apidocs/ee-featurecollection-iterate

[5] MERIT Hydro: Global Hydrography Datasets - Earth Engine Data Catalog https://developers.google.com/earth-engine/datasets/catalog/MERIT_Hydro_v1_0_1

[6] 集水區 - 維基百科中文 https://zh.wikipedia.org/zh-hans/集水區

[7] Nezhad, S. G., Mokhtari, A. R., & Rodsari, P. R. (2017). The true sample catchment basin approach in the analysis of stream sediment geochemical data. Ore Geology Reviews, 83, 127–134. https://doi.org/10.1016/j.oregeorev.2016.12.008

[8] 匯流問題示意圖 http://spatial-analyst.net/ILWIS/htm/ilwisapp/flow_accumulation_functionality.htm

[9] 使用Google Earth Engine計算匯流 - GIS Stack Exchange https://gis.stackexchange.com/questions/279793/calculating-flow-accumulation-using-google-earth-engine

[10] 匯流演算法的GEE分享連結 https://code.earthengine.google.com/ba976e61fe8cc75062494750d99edb10

[11] WWF HydroSHEDS Drainage Direction, 30 Arc-Seconds - Earth Engine Data Catalog https://developers.google.com/earth-engine/datasets/catalog/WWF_HydroSHEDS_30DIR

[12] 批次輸出GEE結果的模板 - GIS Stack Exchange https://gis.stackexchange.com/questions/342461/writing-a-loop-in-google-earth-engine-and-save-all-files-as-unique

相關文章