rgee学习笔记之Image（二）

前面rgee学习笔记之Image（一）介绍了rgee中Image的展示，以及元数据等基本影像信息的获取。下面继续学习Image的用法。

Mathematical operations数值运算

这部分代码主要是进行rgee的波段运算。

• 加载Landsat7数据

# Load two 5-year Landsat 7 composites.
landsat1999 <- ee$Image("LANDSAT/LE7_TOA_5YEAR/1999_2003")
landsat2008 <- ee$Image("LANDSAT/LE7_TOA_5YEAR/2008_2012")

• 使用复杂的方法计算NDVI

# Compute NDVI the hard way.
ndvi1999 <- landsat1999$select("B4")$subtract(landsat1999$select("B3"))$divide(landsat1999$select("B4")$add(landsat1999$select("B3")))

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

# Compute NDVI the easy way.
ndvi2008 <- landsat2008$normalizedDifference(c("B4", "B3"))

• Image.subtract(image2)GEE内置函数，相减

# Compute the multi-band difference image.
diff <- landsat2008$subtract(landsat1999)

• Image.pow(image2)GEE内置函数，幂运算

# Compute the squared difference in each band.
squaredDifference <- diff$pow(2)

• 计算结果的显示

Map$setZoom(zoom = 3)
Map$addLayer(
  eeObject = diff,
  visParams = list(bands = c("B4", "B3", "B2"), min = -32, max = 32),
  name = "difference"
) +
  Map$addLayer(
    eeObject = squaredDifference,
    visParams = list(bands = c("B4", "B3", "B2"), max = 1000),
    name = "squared diff."
  )

波段计算结果

Relational, conditional and Boolean operations关系、条件和布尔运算

• Image.lt(image2)当image1 < image2时，输出1
• 输入是image
• 输出是布尔值
• Image.and(image2)当image1和image2都不为0时，返回1
• 输入是image
• 输出是布尔值

# Load a Landsat 8 image.
image <- ee$Image("LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318")

# Create NDVI and NDWI spectral indices.
ndvi <- image$normalizedDifference(c("B5", "B4"))
ndwi <- image$normalizedDifference(c("B3", "B5"))

# Create a binary layer using logical operations.
bare <- ndvi$lt(0.2)$And(ndwi$lt(0))

# Mask and display the binary layer.
Map$setCenter(lon = -122.3578, lat = 37.7726)
Map$setZoom(zoom = 12)

Map$addLayer(
  eeObject = bare$updateMask(bare),
  visParams = list(min = 0, max = 20),
  name = "bare"
)

黑色部分为运算结果，NDVI<0.2且NDWI<0

Convolutions卷积

ee.Kernel.square

• `ee.Kernel.square(radius, units, normalize, magnitude)
• radius,浮点型，核半径
• units，字符串，默认"pixels"，可以是"pixels"或者"meters"，指定核半径单位，如果核半径单位是米，当缩放级别（zoom-level）变化的时候会变化
• normalize，布尔值，默认："true"，标准化核值和为1（Normalize the kernel values to sum to 1）
• magnitude，浮点型，默认1，确定每个值数量级别（Scale each value by this amount）

ee.Image.convolve

• Image.convolve(kernel) 利用给出的核（Kernel）进行卷积运算
• 图像卷积操作的目的是利用像素点和其邻域像素之前的空间关系，通过加权求和的操作，实现模糊（blurring），锐化（sharpening），边缘检测（edge detection）等功能。
• 卷积的概念和常见运算可以查阅参考文献5

ee.Kernel.laplacian

• `ee.Kernel.laplacian8(magnitude, normalize)``
• 生成一个3×3的Laplacian-8边缘检测核

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

# Define a boxcar or low-pass kernel.
# boxcar = ee$Kernel$square(list(
#   radius = 7, units = 'pixels', 'normalize' = T
# ))

boxcar <- ee$Kernel$square(7, "pixels", T)

# Smooth the image by convolving with the boxcar kernel.
smooth <- image$convolve(boxcar)

# Define a Laplacian, or edge-detection kernel.
laplacian <- ee$Kernel$laplacian8(1, F)

# Apply the edge-detection kernel.
edgy <- image$convolve(laplacian)

Map$setCenter(lon = -121.9785, lat = 37.8694)
Map$setZoom(zoom = 11)

Map$addLayer(
  eeObject = image,
  visParams = list(bands = c("B5", "B4", "B3"), max = 0.5),
  name = "input image"
) +
  Map$addLayer(
    eeObject = smooth,
    visParams = list(bands = c("B5", "B4", "B3"), max = 0.5),
    name = "smoothed"
  ) +
  Map$addLayer(
    eeObject = edgy,
    visParams = list(bands = c("B5", "B4", "B3"), max = 0.5),
    name = "edges"
  )

原始影像

平滑后的影像

锐化（边缘增强）后的影像

Morphological Operations形态操作

ee.Kernel.circle

• `ee.Kernel.circle(radius, units, normalize, magnitude)``
• 生成一个圆形的布尔核

ee.Image.focal_min

• 使用一个定义的核，应用一个形态学的滤波
• `Image.focal_min(radius, kernelType, units, iterations, kernel)``
• radius, 浮点型，默认1.5
• kernelType，字符串，默认“circle”，可选的包括：'circle', 'square', 'cross', 'plus', octagon' and 'diamond'
• units,字符串，默认"pixels"，可选"meters"或"pixels"
• iterations，整型，默认1，应用核的次数
• kernel，核类型，默认：空，一个自定义的核，如果定义了核，kernelType和radius可以忽略
• 类似的函数有：
• ee.Image.focal_median
• ee.Image.focal_mean
• ee.Image.focal_max

# Load a Landsat 8 image, select the NIR band, threshold, display.
# 加载Landsat8影像，选择近红外波段，大于0.2的部分
image <- ee$Image("LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318")$select(4)$gt(0.2)

# Define a kernel.
kernel <- ee$Kernel$circle(radius = 1)

# Perform an erosion followed by a dilation, display.
opened <- image$focal_min(kernel = kernel, iterations = 2)$focal_max(kernel = kernel, iterations = 2)

Map$setCenter(lon = -122.1899, lat = 37.5010)
Map$setZoom(zoom = 13)

Map$addLayer(
  eeObject = image,
  visParams = list(min = 0, max = 1),
  name = "NIR threshold"
) +
  Map$addLayer(
    eeObject = opened,
    visParams = list(min = 0, max = 1),
    name = "opened"
  )

NIR阈值计算结果

核函数形态学滤波后

Gradients梯度

• Image.gradient()计算x,y梯度
• 梯度一般用于增强图像对比度
• 梯度是一个向量
• 梯度的方向是函数变化最快的方向
• 图像边缘处变化大，灰度值变化大，梯度大
• 图像平滑处变化小，灰度值变化小，梯度小

# Load a Landsat 8 image and select the panchromatic band.
image <- ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140318")$select("B8")

# Compute the image gradient in the X and Y directions.
xyGrad <- image$gradient()

# Compute the magnitude of the gradient.
gradient <- xyGrad$select("x")$pow(2)$add(xyGrad$select("y")$pow(2))$sqrt()

# Compute the direction of the gradient.
direction <- xyGrad$select("y")$atan2(xyGrad$select("x"))

# Display the results.
Map$setCenter(lon = -122.054, lat = 37.7295)
Map$setZoom(zoom = 10)

Map$addLayer(
  eeObject = direction,
  visParams = list(min = -2, max = 2),
  name = "direction"
) +
  Map$addLayer(
    eeObject = gradient,
    visParams = list(min = -7, max = 7),
    name = "opened"
  )

梯度方向

梯度数量

参考文献

1. https://csaybar.github.io/rgee-examples/
2. https://developers.google.com/earth-engine/api_docs/#eeimagenormalizeddifference
3. https://developers.google.com/earth-engine/api_docs/#eeimagesubtract
4. https://developers.google.com/earth-engine/api_docs/#eekernelsquare
5. https://blog.csdn.net/u012005313/article/details/84068337
6. https://developers.google.com/earth-engine/api_docs/#eeimagefocal_min
7. https://zhuanlan.zhihu.com/p/113397988