0

我想使用 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'
)
#
4

0 回答 0