在 rgee 中为 COPERNICUS/S2_SR 创建地形校正函数

编程问答 2022-05-01

问题描述

我想使用 NASA SRTM Digital Elevation 30m for DEM 在 S2_SR (BOA) 中创建用于地形校正的 r 函数。我在 GEE 原生函数 (https://mygeoblog.com/2018/07/27/sentinel-2-terrain-correction/) 中找到了一些步骤，但我的想法是将 r 转换为 topoCorr_IC_L4toL8 函数。拜托，我需要一些代码修改的帮助和想法。在我的示例中，使用 BOA 图像去除了云和阴影以及地形校正的完整管道：

在我的例子中：

# Packages -------------------------------------------------------------------
library(raster)
library(rgee)
library(sf)
library(rdgal)

# Initialize the Earth Engine session -----------------------------------------
ee_Initialize(drive = TRUE)

# Selection on interest area 
download.file(
  "https://github.com/Leprechault/trash/raw/main/stands_example.zip",zip_path <- tempfile(fileext = ".zip")
)
unzip(zip_path,exdir = tempdir())

# Selection stands
setwd(tempdir())
roi <- st_read("stands_target.shp") %>% 
  st_bBox() %>% 
  st_as_sfc() %>% 
  sf_as_ee()

# Sentinel-2 MSI dataset into the Earth Engine’s public data archive ------------              
s2 <- ee$ImageCollection("copERNICUS/S2_SR")

# NASA SRTM Digital Elevation 30m -----------------------------------------------

dem <- ee$Image("USGS/SRTMGL1_003")

# Function for topografic correction  --------------------------------------------
topoCorr_IC_L4toL8 <- function(img) {
    # Extract image Metadata about solar position
    var SZ_rad = ee$Image.constant(ee.Number(img.get('MEAN_SOLAR_ZENITH_ANGLE'))).multiply(3.14159265359).divide(180).clip(img.geometry().buffer(10000))
    var SA_rad = ee$Image.constant(ee.Number(img.get('MEAN_SOLAR_AZIMUTH_ANGLE')).multiply(3.14159265359).divide(180)).clip(img.geometry().buffer(10000))
    # Creat terrain layers
    var slp = ee.Terrain.slope(dem).clip(img.geometry().buffer(10000));
    var slp_rad = ee.Terrain.slope(dem).multiply(3.14159265359).divide(180).clip(img.geometry().buffer(10000))
    var asp_rad = ee.Terrain.aspect(dem).multiply(3.14159265359).divide(180).clip(img.geometry().buffer(10000))
    # Calculate the Illumination Condition (IC)
    # Slope part of the illumination condition
    var cosZ = SZ_rad.cos();
    var cosS = slp_rad.cos();
    var slope_illumination = cosS.expression("cosZ * cosS",{'cosZ': cosZ,'cosS': cosS.select('slope')})
    # Aspect part of the illumination condition
    var sinZ = SZ_rad.sin();
    var sinS = slp_rad.sin();
    var cosAziDiff = (SA_rad.subtract(asp_rad)).cos()
    var aspect_illumination = sinZ.expression("sinZ * sinS * cosAziDiff",{'sinZ': sinZ,'sinS': sinS,'cosAziDiff': cosAziDiff});
    # Full illumination condition (IC)
    var ic = slope_illumination.add(aspect_illumination)
    # Add IC to original image
    var img_plus_ic = ee.Image(img.addBands(ic.rename('IC')).addBands(cosZ.rename('cosZ')).addBands(cosS.rename('cosS')).addBands(slp.rename('slope')))
    return (img_plus_ic)
}


# Function for remove cloud and shadows ------------------------------------------
getQABits <- function(image,qa) {
  # Convert decimal (character) to decimal (little endian)
  qa <- sum(2^(which(rev(unlist(strsplit(as.character(qa),"")) == 1))-1))
  # Return a single band image of the extracted QA bits,giving the qa value.
  image$bitwiseAnd(qa)$lt(1)
}

s2_clean <- function(img) {
  # Select only band of interest,for instance,B2,B3,B4,B8
  img_band_selected <- img$select("B[2-12]")
  
  # quality band
  ndvi_qa <- img$select("QA60")

  # Select pixels to mask
  quality_mask <- getQABits(ndvi_qa,"110000000000")
  
  # Mask pixels with value zero.
  img_band_selected$updateMask(quality_mask)
}

# Select S2 images ---------------------------------------------------------------
s2_roi  <- s2$
  filterBounds(roi)$
  filter(ee$Filter$lte("CLOUDY_PIXEL_PERCENTAGE",1))$
  filter(ee$Filter$date('2020-03-01','2021-05-01'))$
  map(s2_clean)

# Get the dates and IDs of the selected images ------------------------------------
nimages <- s2_roi$size()$getInfo()
nimages 
ic_date <- ee_get_date_ic(s2_roi)
ic_date

# Download the results
s2_ic_local <- ee_imagecollection_to_local(
  ic = s2_roi,scale = 10,region = roi,via = 'drive'
)
#

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

rrrgeesatellite-imagesentinel2