rgee学习笔记之rgee get started

rgee作者之一csaybar写的教程，一点点啃……

hello world

下面的代码就是仿照编程最基础的hello world开始

library(rgee)
#ee_reattach() # reattach ee as a reserved word

ee_Initialize()

# traditional R character
print("Hello world!")

# Earth Engine Python Style
ee$String("Hello World from Earth Engine!")$getInfo()
ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140318")$getInfo()

# Earth Engine Pipes Style
ee$String("Hello World from Earth Engine!") %>%
  ee$String$getInfo()

ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140318") %>%
  ee$Image$getInfo()

• ee$String是GEE的字符串，可以自己输入，比如样例中的"hello world"
• ee$Image是GEE的影像对象，可以从GEE服务器中获取
• $getInfo()获取信息

数据加载到地图上

在rgee中，使用leaflet动态地图进行数据的显示交互

加载栅格图

在rgee中要想进行地图交互，需要

1. 导入数据Map$addLayer
2. 定义可视化效果bands/min/max/gamma都要有
3. 把img和vizParams参数一块addLayer，启动显示

# Load an image.
image <- ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140318")

# Display the image.
Map$centerObject(image)
Map$addLayer(image, name = "Landsat 8 original image")

# Define visualization parameters in an object literal.
vizParams <- list(
  bands = c("B5", "B4", "B3"),
  min = 5000, max = 15000, gamma = 1.3
)

# 把img，可视化定义的情况一块输入，地图才可以显示
Map$addLayer(image, vizParams, "Landsat 8 False color")

在viewer中交互显示

加载矢量图

GEE中FeatureCollection是矢量，加载矢量代码示例如下：

# Use Map to add features and feature collections to the map. For example,
counties <- ee$FeatureCollection("TIGER/2016/Counties")

Map$addLayer(
  eeObject = counties,
  visParams = vizParams,
  name = "counties"
)

rgee矢量地图显示

查找影像、矢量并加载

影像筛选加载

sort可以进行排序筛选，例如下面代码中，使用云量筛选前后得到的明显不一样

• filterBounds按范围筛选
• filterDate按时间筛选

# 创建一个ImageCollection
collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1")

#创建一个点
point <- ee$Geometry$Point(-122.262, 37.8719)
#指定起始时间
start <- ee$Date("2014-06-01")
finish <- ee$Date("2014-10-01")

#按位置、时间、云量筛选影像
filteredCollection <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
  filterBounds(point)$
  filterDate(start, finish)$
  sort("CLOUD_COVER", TRUE)

first <- filteredCollection$first()

# Define visualization parameters in an object literal.
vizParams <- list(
  bands = c("B5", "B4", "B3"),
  min = 5000,
  max = 15000,
  gamma = 1.3
)

Map$addLayer(first, vizParams, "Landsat 8 image")

未按云量筛选的影像

按云量筛选的影像

矢量筛选

# Load a feature collection.
featureCollection <- ee$FeatureCollection("TIGER/2016/States")

# Filter the collection.
filteredFC <- featureCollection$filter(ee$Filter$eq("NAME", "California"))

# Display the collection.
Map$addLayer(
  eeObject = filteredFC,
  visParams = list(palette = "red"),
  name = "California"
)

筛选出的加利福尼亚州

波段运算

GEE中内置了大量的函数，可以直接用于波段运算。

• normalizedDifference是GEE中的内置函数，查阅参考文献2中的文档可知：
• 计算公式：(first − second) / (first + second)
• 输入：波段名
• 输出：Image
• Image.subtract(image2)GEE内置函数，相减

# This function gets NDVI from Landsat 5 imagery.
getNDVI <- function(image) {
  return(image$normalizedDifference(c("B4", "B3")))
}

# Load two Landsat 5 images, 20 years apart.
image1 <- ee$Image("LANDSAT/LT05/C01/T1_TOA/LT05_044034_19900604")
image2 <- ee$Image("LANDSAT/LT05/C01/T1_TOA/LT05_044034_20100611")

# Compute NDVI from the scenes.
ndvi1 <- getNDVI(image1)
ndvi2 <- getNDVI(image2)

# Compute the difference in NDVI.
ndviDifference <- ndvi2$subtract(ndvi1)

ndviParams <- list(palette = c(
  "#d73027", "#f46d43", "#fdae61", "#fee08b",
  "#d9ef8b", "#a6d96a", "#66bd63", "#1a9850"
))
ndwiParams <- list(min = -0.5, max = 0.5, palette = c("FF0000", "FFFFFF", "0000FF"))

Map$centerObject(ndvi1)
Map$addLayer(ndvi1, ndviParams, "NDVI 1") +
  Map$addLayer(ndvi2, ndviParams, "NDVI 2") +
  Map$addLayer(ndviDifference, ndwiParams, "NDVI difference")

输出结果有三层：

• NDVI1
• NDVI2
• NDVI1-NDVI2

<<< 左右滑动见更多 >>>

Mapping (what to do instead of a for-loop)

这一节教程我没看明白有啥特别的用途，map看起来好像和terra包里面的app函数有点类似，都是用来给栅格数据集做批量运算的，我理解的是这个意思。如果说的不对还请了解的同学留言指正！

# This function gets NDVI from Landsat 8 imagery.
addNDVI <- function(image) {
  return(image$addBands(image$normalizedDifference(c("B5", "B4"))))
}

# Load the Landsat 8 raw data, filter by location and date.
collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
  filterBounds(ee$Geometry$Point(-122.262, 37.8719))$
  filterDate("2014-06-01", "2014-10-01")

# Map the function over the collection.
ndviCollection <- collection$map(addNDVI)

first <- ndviCollection$first()
print(first$getInfo())

bandNames <- first$bandNames()
print(bandNames$getInfo())

Reducing

Reducer简介

ee.Reducer是GEE中一个具有跨时间、空间、波段、数组等多种数据结构的数据整合函数。

• 可以简单统计整合（minimum, maximum, mean, median, standard deviation等）
• 复杂的一些概要（直方图histogram、线性回归linear regression等）

常见用法

• 跨越时间（over time）imageCollection.reduce()
• 跨越空间image.reduceRegion()、image.reduceNeighborhood()
• 跨波段image.reduce()
• 跨要素属性featureCollection.reduceColumns()，或者以aggregate_开头的FeatureCollection方法

实例

计算了5期Landsat8数据，并取中值

# Load a Landsat 8 collection.
collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
  filterBounds(ee$Geometry$Point(-122.262, 37.8719))$
  filterDate("2014-01-01", "2014-12-31")$
  sort("CLOUD_COVER")

# Compute the median of each pixel for each band of the 5 least cloudy scenes.
median <- collection$limit(5)$reduce(ee$Reducer$median())

# Define visualization parameters in an object literal.
vizParams <- list(
  bands = c("B5_median", "B4_median", "B3_median"),
  min = 5000, max = 15000, gamma = 1.3
)

Map$addLayer(
  eeObject = median,
  visParams = vizParams,
  name = "Median image"
)

Landsat8 Median image

图像统计Image statistics

使用reducer方法进行按区域的计算：

# Load and display a Landsat TOA image.
# 加载Landsat8 TOA数据
image <- ee$Image("LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318")
Map$centerObject(image, zoom = 5)   #地图居中显示
Map$addLayer(
  eeObject = image,
  visParams = list(bands = c("B4", "B3", "B2"), max = 30000),
  name = "Landsat 8"
)

# Create an arbitrary rectangle as a region and display it.
# 创建一个矩形并显示
region <- ee$Geometry$Rectangle(-122.2806, 37.1209, -122.0554, 37.2413)
Map$centerObject(region, zoom=10)
Map$addLayer(
  eeObject = region,
  name = "Region"
)

# Get a dictionary of means in the region.  Keys are bandnames.
# 计算矩形内栅格的均值
mean <- image$reduceRegion(
  reducer = ee$Reducer$mean(),
  geometry = region,
  scale = 30
)

print(mean$getInfo())

各个波段的均值计算结果

掩膜Masking

• ee.Image.mask设置图像的掩膜，掩膜范围内影像值会变成0或接近0
• ee.Image.updateMask 当现有的掩膜范围内不是0得时候，更新掩膜

# This function gets NDVI from Landsat 5 imagery.
getNDVI <- function(image) {
  return(image$normalizedDifference(c("B4", "B3")))
}

# Load two Landsat 5 images, 20 years apart.
image1 <- ee$Image("LANDSAT/LT05/C01/T1_TOA/LT05_044034_19900604")
image2 <- ee$Image("LANDSAT/LT05/C01/T1_TOA/LT05_044034_20100611")

# Compute NDVI from the scenes.
ndvi1 <- getNDVI(image1)
ndvi2 <- getNDVI(image2)

# Compute the difference in NDVI.
ndviDifference <- ndvi2$subtract(ndvi1)
# Load the land mask from the SRTM DEM.陆地范围掩膜
landMask <- ee$Image("CGIAR/SRTM90_V4")$mask()

# Update the NDVI difference mask with the land mask.使用陆地范围掩膜NDVI
maskedDifference <- ndviDifference$updateMask(landMask)

# Display the masked result.
vizParams <- list(
  min = -0.5,
  max = 0.5,
  palette = c("FF0000", "FFFFFF", "0000FF")
)

Map$addLayer(
  eeObject = maskedDifference,
  visParams = vizParams,
  name = "NDVI difference"
)

掩膜后的NDVI差值

综合练习

最后是一个综合练习，综合了前面用到的知识点：

• NDVI波段计算normalizedDifference
• 掩膜image$updateMask(clouds$lt(10))
• 应用计算map
• 数据整合Reducer
• 区域统计和结果输出

# This function gets NDVI from Landsat 8 imagery.
addNDVI <- function(image) {
  return(image$addBands(image$normalizedDifference(c("B5", "B4"))))
}

# This function masks cloudy pixels.
cloudMask <- function(image) {
  clouds <- ee$Algorithms$Landsat$simpleCloudScore(image)$select("cloud")
  return(image$updateMask(clouds$lt(10)))
}

# Load a Landsat collection, map the NDVI and cloud masking functions over it.
collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1_TOA")$
  filterBounds(ee$Geometry$Point(c(-122.262, 37.8719)))$
  filterDate("2014-03-01", "2014-05-31")$
  map(addNDVI)$
  map(cloudMask)

# Reduce the collection to the mean of each pixel and display.
meanImage <- collection$reduce(ee$Reducer$mean())
vizParams <- list(
  bands = c("B5_mean", "B4_mean", "B3_mean"),
  min = 0,
  max = 0.5
)

Map$addLayer(
  eeObject = meanImage,
  visParams = vizParams,
  name = "mean"
)

# Load a region in which to compute the mean and display it.
counties <- ee$FeatureCollection("TIGER/2016/Counties")
santaClara <- ee$Feature(counties$filter(
  ee$Filter$eq("NAME", "Santa Clara")
)$first())

Map$addLayer(
  eeObject = santaClara,
  visParams = list(palette = "yellow"),
  name = "Santa Clara"
)

# Get the mean of NDVI in the region.
mean <- meanImage$select("nd_mean")$reduceRegion(
  reducer = ee$Reducer$mean(),
  geometry = santaClara$geometry(),
  scale = 30
)

# Print mean NDVI for the region.
cat("Santa Clara spring mean NDVI:", mean$get("nd_mean")$getInfo())

<<< 左右滑动见更多 >>>

参考文献

1. https://csaybar.github.io/rgee-examples/
2. https://developers.google.com/earth-engine/api_docs/
3. https://developers.google.com/earth-engine/guides/reducers_intro
4. https://developers.google.com/earth-engine/api_docs/#eeimagemask

更多rgee资料请阅读#rgee