SamplePredictorLayers.R

SamplePredictorLayers.R

#' Script to sample multiple layers by years 1990 - 2011 using a point dataset
#' @description Layers: Landsat spectral values; NDVI; NDMI; NBR; Elevation, Slope, Aspect.
#' 
#' PARAMETERS
#' 
#' 
library(rgee)
ee_Initialize()

USUptsNew <- ee$FeatureCollection("users/dougramsey/ForCaMF/PointData/USU_RandomPoints_NEW")
DEM <- ee$Image("USGS/NED")

Year <- 2010
Years <- ee$List$sequence(1990,2011)
Years2 <- c(1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,
            2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,
            2010,2011)
startMonth <- 7
endMonth <- 8

FCPalette <- c("3d4551","ff0000","ffff00","00ffff","ffffff","ff00ff")
FCpixValues <- 1:6
FCclassNames <- c('Undisturbed','Fire','Harvest','Insects','??','Abiotic')

LCPalette <- c("3d4551","f39268","d54309","00a398","1B1716")
LCpixValues <- 1:5
LCclassNames <- c('Stable','Slow Loss','Fast Loss','Gain','Not Processed')

points <- USUptsNew
# region <- '04'
#  restrict the geography.
fc <- points



# ----------------
#  FUNCTION BLOCK-
# ----------------

# Filter for clouds L457
masksr <- function(image) {
  qa <- image$select("pixel_qa")
  # If the cloud bit (5) is set and the cloud confidence (7) is high
  # or the cloud shadow bit is set (3), then it's a cloudy pixel.
  cloud <- qa$bitwiseAnd(2**5)$And(qa$bitwiseAnd(2**7))$Or(qa$bitwiseAnd(2**3))
  # Remove edge pixels that don't occur in all bands
  mask2 <- image$mask()$reduce(ee$Reducer$min())
  image$updateMask(cloud$Not())$updateMask(mask2)
}

# NDVI Functions for L8 and L7 (the L7 function is used for L4&5)
addIndices <- function(image) {
  ndvi <- image$normalizedDifference(c('NIR', 'R'))$rename('NDVI')
  ndmi <- image$normalizedDifference(c('NIR', 'SWIR1'))$rename('NDMI')
  nbr <- image$normalizedDifference(c('SWIR2', 'NIR'))$rename('NBR')
  image$addBands(ndvi)$addBands(ndmi)$addBands(nbr)
}

# ----------------
#  INPUT DATA ----
# ----------------
# This sets up the selection and renaming of bands so 
# all images will have the same spectral bands

l8bands <- c('B2','B3','B4','B5','B6','B7')
l8bndRename <- c('B','G','R','NIR','SWIR1','SWIR2')
l57bands <- c('B1','B2','B3','B4','B5','B7')
l57bndRename <- c('B','G','R','NIR','SWIR1','SWIR2')

# Get the Image Collections and merge into one
# Load the Landsat 8 ImageCollection.  Mask clouds and add NDVI.
l8 <- ee$ImageCollection('LANDSAT/LC08/C01/T1_SR')$
  filterBounds(fc)$
  map(masksr)$
  map(function(image){image$select(l8bands)$rename(l8bndRename)})$
  map(addIndices)

# Load the Landsat 7 ImageCollection.  Mask clouds and add NDVI. 
l7 <- ee$ImageCollection('LANDSAT/LE07/C01/T1_SR')$
  filterBounds(fc)$
  map(masksr)$
  map(function (image){image$select(l57bands)$rename(l57bndRename)})$
  map(addIndices)

# Load the Landsat 5 ImageCollection.  Mask clouds and add NDVI.
l5 <- ee$ImageCollection('LANDSAT/LT05/C01/T1_SR')$
  filterBounds(fc)$
  map(masksr)$
  map(function (image){image$select(l57bands)$rename(l57bndRename)})$
  map(addIndices)

Merged <- l5$merge(l7)$sort("system:time_start")
Merged <- Merged$merge(l8)$sort("system:time_start")$filter(
  ee$Filter$calendarRange(Years$get(0), Years$get(-1), 'year')
)

# print(Merged)

Aspect <- ee$Terrain$aspect(DEM)
Slope <- ee$Terrain$slope(DEM)


# -------------------------------------------------------------------------
#          GENERATE ForCaMF IMAGE COLLECTION FOR ALL YEARS
# Sort the ForCaMF Codes table and extract the DistType and Value attributes 
# for all records and place them in their own lists (set up for remapping the 
# image) Remap the ForCaMF_Mosaic image into separate images for each year 
# -------------------------------------------------------------------------
ForCaMFcodes <- ee$FeatureCollection("users/dougramsey/ForCaMF/MOSA10_combinations")
ForCaMF <- ee$Image("users/dougramsey/ForCaMF/ForCaMF_Mosaic")

#  Sort the codes based on the grid code value
comb <- ForCaMFcodes$sort('Value')

#  Extract the disturbance codes and place them in a list
distType <- comb %>% 
  ee$FeatureCollection$toList(
    comb$size()
  ) %>% 
  ee$List$map(
    ee_utils_pyfunc(
      function(feat){ee$String(ee$Feature(feat)$get('DistType'))}
    )
  )
  
# Extract the grid codes into a separate list which should match, in 
# order the disturbance codes
grdCode <- comb %>% 
  ee$FeatureCollection$toList(comb$size()) %>% 
  ee$List$map(
    ee_utils_pyfunc(
      function(feat){ee$Number$parse(ee$String(ee$Feature(feat)$get('Value')))}      
    )
  )

#  Generate an image collection to hold ForCaMF disturbance images for each year.
ForCaMF_Dist <- ee$ImageCollection$fromImages(
  Years$map(
    ee_utils_pyfunc(
      function (Y) {
        # Extract the disturbance code for the specified year
        yearDist <- distType$map(
          ee_utils_pyfunc(
            function(d) {
              ee$Number$parse(ee$String(d)$split('')$remove('d')$get(Years$indexOf(Y)))
            }
          )
        )
        #  Remap the Mosa_Region10f image to correspond to a specific years disturbance code.
        Dist <- ForCaMF$select('b1')$remap(grdCode, yearDist)
        # New band Name
        y <- ee$String('ForCaMF')
        # return the image, set the year and band name
        Dist$set('year', Y)$rename(y)
      }      
    )
  )
)

# print ('ForCaMF_Dist',ForCaMF_Dist)

# =========================================================================
# =========================================================================
#  EXTRACT LCMS FOR SAME YEARS AS ForCaMF and Mask LCMS by ForCaMF
#  Mask LCMS with ForCaMF and change the band name to 'LCMS_'+year
Mask <- function(image) {
  #  var year = image.get('year');
  image$updateMask(ForCaMF)
}

LCMS <- ee$ImageCollection("USFS/GTAC/LCMS/v2020-5")
LCMS_Dist <- LCMS$filter(ee$Filter$calendarRange(Years$get(0),Years$get(-1),'year'))$
  filter(ee$Filter$eq('study_area', 'CONUS'))$
  select(c('Change','Land_Cover','Change_Raw_Probability_Slow_Loss',
         'Change_Raw_Probability_Fast_Loss','Change_Raw_Probability_Gain'))$map(Mask)
# print('LCMS',LCMS_Dist)
slopeMask <- Slope$updateMask(ForCaMF)
aspectMask <- Aspect$updateMask(ForCaMF)
demMask <- DEM$updateMask(ForCaMF)

# =========================================================================
# =========================================================================
#                      COMBINE LCMS AND ForCaMF
#  Generate multi-band images for each year with first band as ForCaMF 
#  and second band as LCMS, followed by DEM, slope, aspect, and a year band
combined_Dist <- ee$ImageCollection$fromImages(
  Years$map(function (Y) {
    FCimage <- ForCaMF_Dist$filter(ee$Filter$eq('year', Y))$first()
    LCimage <- LCMS_Dist$filter(ee$Filter$eq('year', Y))$first()
    
    theYear <- ee$Image$constant(Y)$uint16()$rename('year')
    theYear <- theYear$updateMask(FCimage$mask())
    
    #  ASSUMPTION ON THE MONTH RANGE.  CHECK ****************************************    
    Landsat <- Merged$filter(ee$Filter$calendarRange(Y, Y, 'year'))$
      filter(ee$Filter$calendarRange(startMonth, endMonth, 'month'))$
      median()$updateMask(FCimage)
      
    FCimage$addBands(LCimage)$addBands(DEM)$
      addBands(slopeMask)$addBands(aspectMask)$
      addBands(theYear)$addBands(Landsat)$set('year',Y)
  }))

# print ('combined_Dist',combined_Dist)

# =========================================================================
# =========================================================================
#  Sample LCMS by stratifying the samples acrosss ForCaMF classses
# =========================================================================
samples <- Years2$map(
  function (Y) {
  FCband <- ee$String('ForCaMF')
  LCband <- ee$String('Change')
  image <- combined_Dist$filterMetadata('year','equals',Y)$first()
  masked <- image
  masked$sampleRegions(
    collection = points,
    scale = 30, 
    projection = 'EPSG:5070', 
    tileScale = 1, 
    geometries = TRUE 
  )}
)

allSamples <- ee$FeatureCollection(samples)$flatten()