antony

rgee.R

##  Travel time to school facilities as a marker of geographical      ## 
##         accessibility across heterogeneous land cov                ##
##' @GabrielCarrasco-Escobar, @EdgarManrique, @KellyTello-Lizarraga,  ##
##'              @J.JaimeMiranda,@AntonyBarja                         ##
##'-------------------------------------------------------------------##

# Requerements
library(tidyverse)
library(cptcity)
library(rgee)
library(sf)

ee_Initialize()

# Preparing dataset (pq cargas una geometria tan pesada!, siempre usa un bbox!!)
peru  <- ee$FeatureCollection('users/edgarmanrique30/Peru_geometry/Limite_departamental') %>% 
  ee$FeatureCollection$geometry() %>% 
  ee$Geometry$bounds() %>% 
  ee_as_sf() %>% 
  st_bbox() %>% 
  ee$Geometry$Rectangle()

# DEM y slope
dem   <- ee$Image('USGS/SRTMGL1_003') %>% 
  ee$Image$clip(peru)
slope <- ee$Terrain$slope(dem)

# Land use
landc <- ee$ImageCollection('MODIS/006/MCD12Q1') %>% 
  ee$ImageCollection$select("LC_Type1") %>% 
  ee$ImageCollection$filterDate("2017-01-01", "2017-12-31") %>% 
  ee$ImageCollection$median() %>% 
  ee$Image$clip(peru)

# vias
vias_dep <- ee$FeatureCollection('users/edgarmanrique30/Peru_geometry/red_vial_departamental_dic16')
vias_nac <- ee$FeatureCollection('users/edgarmanrique30/Peru_geometry/red_vial_nacional_dic16')
vias_vec <- ee$FeatureCollection('users/edgarmanrique30/Peru_geometry/red_vial_vecinal_dic16')

rios  <- ee$Image('WWF/HydroSHEDS/15ACC') %>% 
  ee$Image$gt(5000) %>% 
  ee$Image$remap(c(0,1), c(0,9))

anp   <- ee$FeatureCollection('users/edgarmanrique30/Peru_geometry/ANP-Nacional')
school <- ee$FeatureCollection('users/antonybarja8/school')


black <- ee$Image(0)$byte()
vias_nac <- black$paint(vias_nac, 80)$clip(peru)
vias_dep <- black$paint(vias_dep, 50)$clip(peru)
vias_vec <- black$paint(vias_vec, 30)$clip(peru)

# LC_Type1, Remapping the pixel values of each category of land cover to their respective speed in km/h.
landcspeed <- landc$
  remap(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 ,14, 15, 16, 17),
        c(3.24, 1.62, 3.24, 4, 3.24, 3, 4.2, 4.86, 4.86, 4.86, 2, 2.5, 5, 3.24, 1.62, 3, 1))
# Filtering urban areas to multiply the roads speed by .7
landc_urban <- landc$eq(13) 
# Filtering urban areas to multiply the roads speed by .7
landc_urban <- landc_urban$remap(c(0,1),c(1,0.7)) 
vias_nac <- vias_nac$multiply(landc_urban) # Multiplying roads layers by 0.7 in urban areas
vias_dep <- vias_dep$multiply(landc_urban) # Multiplying roads layers by 0.7 in urban areas
vias_vec <- vias_vec$multiply(landc_urban) # Multiplying roads layers by 0.7 in urban areas

# Creating Natural protected areas layer
anp <- black$paint(anp, 1)$
  clip(peru) 
#Remapping values to 0.2 km/h of Natural protected areas to multiply Landcover speed
anp <- anp$remap(c(0,1),c(1,0.2)) 
landcspeed <- landcspeed$multiply(anp) # Multiplying Landcover speed by 0.2 on Natural protected areas
landcspeed <- landcspeed$toDouble()$select(list(0),list("speed"))
rios <- rios$toDouble()$select(list(0),list("speed"))
vias_nac <- vias_nac$toDouble()$select(list(0),list("speed")) # unifying the band name
vias_dep <- vias_dep$toDouble()$select(list(0),list("speed")) # unifying the band name
vias_vec <- vias_vec$toDouble()$select(list(0),list("speed")) # unifying the band name
# Mergging all layers into a collection
collection <- ee$ImageCollection(
  list(landcspeed,
       rios,
       vias_nac,
       vias_dep,
       vias_vec
  )
) 

fsurface <- collection$max() # Calculating the maximum value of speed on a single pixel
# eaf <- function(x) {1.01*exp(-0.0001072*x)} # Elevation adjustment factor
eaf <- dem$expression(
  c('1.01*exp(-0.0001072*DEM)'),list(
    'DEM'= dem$select('elevation')))


# thf <- function(x) {6*exp(-3.5*abs((tan(x/57.296) + 0.05)))/5} # Tobler's hikking function adjustment
thf <- slope$expression(
  c('6*exp(-3.5*abs((tan(slope/57.296) + 0.05)))/5'), list(
    'slope'= slope$select(list(0))
  ))

# Adjusting the friction surface by EAF and THF
fsurface <- fsurface$multiply(eaf)$multiply(thf)

# convert <- function(x) {(x * 1000 / 60) ^ -1} # converts km/h to min/m
fsurface <- fsurface$expression(
  c('(x * 1000 / 60) ** -1'),
  list('x'= fsurface$select(list(0)))
)

# Paint the input points, essentially converting them to a raster.
# Theoretically this will merge any points that fall within the same pixel (of the resulting 30-arc-second resolution).
sources <- black$paint(school, 1)
sources <- sources$updateMask(sources)
# Compute the min cost path distance map, with a horizon of 1500 km.
# This can be reduced for high-latitude areas and/or to shorten processing time.
distance <- fsurface$cumulativeCost(sources, 400000)  # The function accepts meters rather than km.
distance <- ee$Image(distance)$toInt() # Here we convert the output to integer to make the output .tif smaller (and the code more likely to run successfully).
distance <- distance$clip(peru)$reproject("EPSG:4326")

# Final Map

viz <- list(
  min = 0,
  max = 600,
  palette = cpt(pal = 'grass_bcyr',n = 5,rev = T)
)

Map$addLayer(distance,viz,'accessibility') + 
  Map$addLayer(fsurface$clip(peru),viz,'friction') + 
  Map$addLayer(school)

# Extract value of IIEE
extract_value <- function(x, y, by = 1000){
  y_len <- y$size()$getInfo()
  for(i in seq(1,y_len,by)) { 
    index <-  i - 1
    print(sprintf("Extracting information [%s/%s]...",index,y_len))
    ee_value_school <- ee$FeatureCollection(y) %>% 
      ee$FeatureCollection$toList(by, index) %>% 
      ee$FeatureCollection()
    
  if(i == 1) {
      dataset <- ee_extract(
        x = x,
        fun = ee$Reducer$median(),
        y = ee_value_school,
        scale = 10000, 
        sf = F)
      } else {
      db_local <- ee_extract(x = x,y = ee_value_school,
                             fun = ee$Reducer$median(),
                             sf = TRUE)
      dataset <- rbind(da(taset,db_local)
    }
  }
  return(dataset)
}  

school_value <- extract_value(x = distance, y = school)
school %>% ee_get(0) %>% ee$FeatureCollection$getInfo()
Map$addLayer(eeObject = distance, visParams = list(min = 0, max = 1))
rr1 = ee_as_raster(distance, scale = 10000)