Начнем с вычисления, необходимого для создания полосы, которая показывает пиксели, где данные Хансена и др. показывают как потерю, так и прирост.
В наборе данных Хансена и др. есть полоса, пиксели которой равны 1, если произошла потеря, и 0 в противном случае (loss), и полоса, которая равна 1, если произошел прирост, и 0 в противном случае (gain). Чтобы создать полосу, в которой пиксели в полосах loss и gain имеют значение 1, вы можете использовать логический метод and() для изображений. Метод and() вызывается как image1.and(image2) и возвращает изображение, в котором пиксели равны 1, если и image1, и image2 равны 1, и 0 в остальных случаях:
Редактор кода (JavaScript)
// Load the data and select the bands of interest. var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015'); var lossImage = gfc2014.select(['loss']); var gainImage = gfc2014.select(['gain']); // Use the and() method to create the lossAndGain image. var gainAndLoss = gainImage.and(lossImage); // Show the loss and gain image. Map.addLayer(gainAndLoss.updateMask(gainAndLoss), {palette: 'FF00FF'}, 'Gain and Loss');
Результат, увеличенный до Арканзаса в режиме спутника, должен выглядеть примерно так, как показано на рисунке 1.
Объединив этот пример с результатом из предыдущего раздела, можно воссоздать рисунок из начала руководства:
Редактор кода (JavaScript)
// Displaying forest, loss, gain, and pixels where both loss and gain occur. var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015'); var lossImage = gfc2014.select(['loss']); var gainImage = gfc2014.select(['gain']); var treeCover = gfc2014.select(['treecover2000']); // Use the and() method to create the lossAndGain image. var gainAndLoss = gainImage.and(lossImage); // Add the tree cover layer in green. Map.addLayer(treeCover.updateMask(treeCover), {palette: ['000000', '00FF00'], max: 100}, 'Forest Cover'); // Add the loss layer in red. Map.addLayer(lossImage.updateMask(lossImage), {palette: ['FF0000']}, 'Loss'); // Add the gain layer in blue. Map.addLayer(gainImage.updateMask(gainImage), {palette: ['0000FF']}, 'Gain'); // Show the loss and gain image. Map.addLayer(gainAndLoss.updateMask(gainAndLoss), {palette: 'FF00FF'}, 'Gain and Loss');
Количественная оценка изменений лесного покрова в интересующем регионе
Теперь, когда вы лучше знакомы с каналами в наборе данных Хансена и др., мы можем использовать изученные концепции для расчета статистики прироста и потери лесов в интересующем регионе. Для этого нам понадобятся векторные данные (точки, линии и многоугольники). Векторный набор данных представлен в Earth Engine как FeatureCollection.
Подробнее о коллекциях объектов и о том, как импортировать векторные данные…
В этом разделе мы сравним общую площадь лесов, вырубленных в Республике Конго в 2012 году, с площадью лесов, вырубленных на охраняемых территориях страны за тот же период.
Как вы узнали из руководства по Earth Engine API, ключевой метод для расчета статистики в регионе изображения – reduceRegion(). Подробнее о том, как уменьшить количество регионов… Например, предположим, что мы хотим рассчитать количество пикселей, которые, по оценкам, представляют собой потерю леса за период исследования. Для этого используйте следующий код:
Редактор кода (JavaScript)
// Load country features from Large Scale International Boundary (LSIB) dataset. var countries = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017'); // Subset the Congo Republic feature from countries. var congo = countries.filter(ee.Filter.eq('country_na', 'Rep of the Congo')); // Get the forest loss image. var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015'); var lossImage = gfc2014.select(['loss']); // Sum the values of forest loss pixels in the Congo Republic. var stats = lossImage.reduceRegion({ reducer: ee.Reducer.sum(), geometry: congo, scale: 30 }); print(stats);
В этом примере используется редуктор ee.Reducer.sum(), чтобы суммировать значения пикселей в lossImage в пределах объекта congo. Поскольку lossImage состоит из пикселей, имеющих значение 1 или 0 (для потери или отсутствия потери соответственно), сумма этих значений эквивалентна количеству пикселей потери в регионе.
Если запустить скрипт в таком виде, возникнет ошибка:
По умолчанию максимальное количество пикселей в reduceRegion() составляет 10 миллионов. Это сообщение об ошибке указывает на то, что Республика Конго занимает около 383 миллионов пикселей Landsat. К счастью, у функции reduceRegion() есть много параметров, один из которых (maxPixels) позволяет контролировать, сколько пикселей используется в вычислениях. Если указать этот параметр, вычисление будет выполнено успешно:
Редактор кода (JavaScript)
// Load country features from Large Scale International Boundary (LSIB) dataset. var countries = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017'); // Subset the Congo Republic feature from countries. var congo = countries.filter(ee.Filter.eq('country_na', 'Rep of the Congo')); // Get the forest loss image. var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015'); var lossImage = gfc2014.select(['loss']); // Sum the values of forest loss pixels in the Congo Republic. var stats = lossImage.reduceRegion({ reducer: ee.Reducer.sum(), geometry: congo, scale: 30, maxPixels: 1e9 }); print(stats);
Разверните объект, выведенный в консоль, и посмотрите результат. Он должен быть равен 4 897 933 пикселям. Вы можете немного очистить распечатку в консоли, пометив выходные данные и получив интересующий результат из словаря, возвращенного reduceRegion():
Редактор кода (JavaScript)
print('pixels representing loss: ', stats.get('loss'));
Расчет площади пикселей
Вы почти готовы ответить на вопрос о том, какая площадь была потеряна в Республике Конго и какая часть из нее находилась на охраняемых территориях. Осталось преобразовать пиксели в реальную площадь. Это преобразование важно, поскольку мы не всегда знаем размер пикселей, передаваемых в reduceRegion(). Для вычисления площади в Earth Engine есть метод ee.Image.pixelArea(), который генерирует изображение, где значение каждого пикселя – это его площадь в квадратных метрах. Умножив изображение потерь на изображение площади и сложив результаты, мы получим меру площади:
Редактор кода (JavaScript)
// Load country features from Large Scale International Boundary (LSIB) dataset. var countries = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017'); // Subset the Congo Republic feature from countries. var congo = countries.filter(ee.Filter.eq('country_na', 'Rep of the Congo')); // Get the forest loss image. var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015'); var lossImage = gfc2014.select(['loss']); var areaImage = lossImage.multiply(ee.Image.pixelArea()); // Sum the values of forest loss pixels in the Congo Republic. var stats = areaImage.reduceRegion({ reducer: ee.Reducer.sum(), geometry: congo, scale: 30, maxPixels: 1e9 }); print('pixels representing loss: ', stats.get('loss'), 'square meters');
Теперь результат составляет 4 372 575 052 квадратных метров, потерянных за период исследования.
Теперь вы готовы ответить на вопрос, заданный в начале этого раздела: какая площадь лесов была утрачена в Республике Конго в 2012 году и какая часть этих лесов находилась на охраняемых территориях?
Редактор кода (JavaScript)
// Load country features from Large Scale International Boundary (LSIB) dataset. var countries = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017'); // Subset the Congo Republic feature from countries. var congo = ee.Feature( countries .filter(ee.Filter.eq('country_na', 'Rep of the Congo')) .first() ); // Subset protected areas to the bounds of the congo feature // and other criteria. Clip to the intersection with congo. var protectedAreas = ee.FeatureCollection('WCMC/WDPA/current/polygons') .filter(ee.Filter.and( ee.Filter.bounds(congo.geometry()), ee.Filter.neq('IUCN_CAT', 'VI'), ee.Filter.neq('STATUS', 'proposed'), ee.Filter.lt('STATUS_YR', 2010) )) .map(function(feat){ return congo.intersection(feat); }); // Get the loss image. var gfc2014 = ee.Image('UMD/hansen/global_forest_change_2015'); var lossIn2012 = gfc2014.select(['lossyear']).eq(12); var areaImage = lossIn2012.multiply(ee.Image.pixelArea()); // Calculate the area of loss pixels in the Congo Republic. var stats = areaImage.reduceRegion({ reducer: ee.Reducer.sum(), geometry: congo.geometry(), scale: 30, maxPixels: 1e9 }); print( 'Area lost in the Congo Republic:', stats.get('lossyear'), 'square meters' ); // Calculate the area of loss pixels in the protected areas. var stats = areaImage.reduceRegion({ reducer: ee.Reducer.sum(), geometry: protectedAreas.geometry(), scale: 30, maxPixels: 1e9 }); print( 'Area lost in protected areas:', stats.get('lossyear'), 'square meters' );
Результат показывает, что из 348 036 295 квадратных метров леса, потерянных в Республике Конго в 2012 году, 11 880 976 квадратных метров находились на охраняемых территориях, как показано в таблице Всемирной базы данных по охраняемым территориям.
Единственные изменения между этим скриптом и предыдущим – это добавление информации об охраняемой территории и изменение скрипта с просмотра общих потерь на просмотр потерь в 2012 году. Для этого потребовалось два изменения. Во-первых, это новое изображение lossIn2012, в котором пиксели со значением 1 соответствуют местам, где в 2012 году была зафиксирована потеря лесного покрова, а пиксели со значением 0 – местам, где потери не было. Во-вторых, поскольку название группы отличается (lossyear вместо loss), название свойства в операторе print пришлось изменить.
В следующем разделе мы рассмотрим некоторые продвинутые методы расчета и построения диаграмм потери лесов за каждый год, а не только за один год, как мы это делали в этом разделе.