cloud-sorting and visualizing satellite data

setup:

library(sf)
library(rgee)
library(dplyr)
library(ggplot2)
library(plotly)
library(mapview)

ee_Initialize()

## -- rgee 1.0.5 ------ earthengine-api 0.1.233 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/kudoh 
## -----------------------------------------------

mapviewOptions(basemaps = c("CartoDB.DarkMatter", "Esri.WorldImagery"))

study points:

point = st_point(c(-84.35128880352657,46.491522302001066)) %>% sf_as_ee()

geom = st_point(c(-86.12798881609389,45.56853951395261)) %>% sf_as_ee()

geom2 = st_point(c(-85.20019149611451,45.98488673090658)) %>% sf_as_ee()

geom3 = st_point(c(-85.02672767470338,45.545523529351975)) %>% sf_as_ee()

geom4 = st_point(c(-86.62875928018774,45.72010451925803)) %>% sf_as_ee()

part i: scanning for cloud-free images

scene = ee$ImageCollection('LANDSAT/LC08/C01/T1_TOA')$ # getting landsat imagery 
  filterBounds(point)$ # spatial filter
  filterDate('2015-01-01', '2015-12-31')$ # temporal filter
  sort('CLOUD_COVER') %>%  # sorting by cloud cover
  {.$toList(.$size())} %>% # making a list of images
  {ee$Image(.$get(0))} # getting first image from list

Map$centerObject(scene, 9)

part ii: true-color, false color, and natural color composites

truecolor = list(
  bands = c('B4', 'B3', 'B2'),
  min = c(
    -0.14659746099527587,-0.14659746099527587,-0.14659746099527587
  ),
  max = c(
    0.40796423639633567,
    0.40796423639633567,
    0.40796423639633567
  )
)

falsecolor = list(
  bands = c('B5', 'B4', 'B3'),
  min = c(
    -0.17543534027220278,
    -0.17543534027220278,
    -0.17543534027220278
  ),
  max = c(
    0.5290458963180128,
    0.5290458963180128,
    0.5290458963180128
  )
)

naturalcolor = list(
  bands = c('B7', 'B5', 'B3'),
  min = c(
    -0.1437304446502991,
    -0.1437304446502991,
    -0.1437304446502991
  ),
  max = c(
    0.3671647239998674,
    0.3671647239998674,
    0.3671647239998674
  )
)

up = ee_as_mapview(
  Map$addLayer(scene, truecolor, 'true-color composite') +
    Map$addLayer(scene, falsecolor, 'false-color composite') +
    Map$addLayer(scene, naturalcolor, 'natural-color composite')
)

part iii: spectral characteristics of features

scene = scene$select('B2', 'B3', 'B4', 'B5', 'B6', 'B7')

point1 = ee$Feature(point, list(label = "town"))
point2 = ee$Feature(geom, list(label = "water"))
point3 = ee$Feature(geom2, list(label = "ice/snow"))
point4 = ee$Feature(geom3, list(label = "woods"))
point5 = ee$Feature(geom4, list(label = "farmland"))

points = ee$FeatureCollection(list(point1,point2,point3,point4,point5))

wv = c(.48, .54, .65, .86, 1.6, 2.2) # mean wavelengths of bands

# either $sampleRegions() or ee$Reducer$first() seem to work, though ee$Reducer$first() was used in the original lab 

# bv = scene$sampleRegions( # band values
#   collection = points,
#   properties = list('label'),
#   scale =30
# )$getInfo()$features

bv = scene$reduceRegions( # band values
  collection = points,
  reducer = ee$Reducer$first(),
  scale =1
)$getInfo()$features

spectral = tibble(band = character(),
                  feature = character(),
                  `brightness value` = double())

for (i in seq_along(bv)){
 for (j in 1:6){
   row = list(
   band = bv[[i]][['properties']][j] %>%  names(),
   feature  =  bv[[i]][['properties']][['label']],
   "brightness value" = bv[[i]][['properties']][[j]] %>%  round(3)
   )
   spectral = bind_rows(spectral, row)
 }
}

spectral = mutate(
  spectral,
  wavelength = case_when(
    band == 'B2' ~ wv[1],
    band == 'B3' ~ wv[2],
    band == 'B4' ~ wv[3],
    band == 'B5' ~ wv[4],
    band == 'B6' ~ wv[5],
    band == 'B7' ~ wv[6]
  )
)
  
sv = ggplotly(spectral %>% 
  ggplot(aes(x = wavelength, y =  `brightness value`, color = feature)) + 
  geom_line()+
  geom_point() + 
  labs(title = 'Spectra values in sample points',
      y = 'Brightness Values',
      x = 'Bands') +
  theme(
    legend.title = element_blank()
  )
)