## fix sampling points / location from sf data
th_sf_sample <-
  th_sf %>%
  group_by(NAME_1) %>%
  nest() %>%
  # Make grids on the map by province
  dplyr::mutate(
    province_grid_centers = furrr::future_map(
      data,
      ~
        # Make grids
        sf::st_make_grid(
          x = .,
          cellsize = cellsize,
          what = "centers"
        ) %>%
        # Obtain crossing point of the grids
        sf::st_intersection(.) %>%
        # Bind the coordinates of crossing points.
        base::do.call(rbind, .) %>%
        dplyr::as_tibble() %>%
        stats::setNames(c("lon","lat")
        )
    )
  ) %>%
  # Sample the coordinates of crossing points by province
  dplyr::mutate(
    sample_grid = furrr::future_map2(
      .x = province_grid_centers,
      .y = sample_size,
      .f =
        ~
        dplyr::sample_n(.x,
                        size = .y,
                        replace = FALSE
        )
    )
  ) %>%
  tidyr::unnest(sample_grid) %>%
  dplyr::bind_cols(
    individual = factor(
      rep(
        c(
          1:sample_size
        ),
        times = length(th_sf$NAME_1)
      )
    )
  )
# ## save all the results
# Note
# Comment out when not in use. The procedure below takes long computation period.
# # This process needs extremely long periods.
# saveRDS(th_sf_sample, "th_sf_sample.rds")
# gc()
# gc()
# gc()
#

# Obtain satellite imagery from Google maps
#

th_sf_sample <- readRDS("th_sf_sample.rds")
#
th_sf_rgb <-
  # th_sf_sample$sample_grid %>%
  # # bind the grids' coordinates together
  # base::do.call(dplyr::bind_rows,.) %>%
  # obtain the satellite imagery
  th_sf_sample %>%
  dplyr::mutate(
    ggmap_image_satellite = purrr::map2(
      lon,
      lat,
      function(lon, lat)
      {
        ggmap::get_map(
          location = c(
            lon = lon,
            lat = lat
          ),
          maptype = "satellite",
          size = c(320, 320),
          scale = 1,
          format = "png",
          zoom = 18 # <- BEWARE!!
        )
      }
    )
  ) %>%
  # Note on factor number and aggregate function.
  # The factor numbers depend on our previous experiments.
  # 2, 4, and 8 need too long computation period 
  dplyr::mutate(
    fact_number_10 = c(10), # dense
    fact_number_20 = c(20),  # modest
    fact_number_20 = c(40)  # rough
  ) %>%
  tidyr::gather(key = "fn",  value = "fact_number", starts_with("fact_number_")) %>%
  # Note on the aggregate function
  # The function, "mean", depends on our previous experiments.
  mutate(
    aggregate_fun_mean = c("mean")
  ) %>%
  tidyr::gather(
    key = "af",  
    value = "aggregate_fun", 
    starts_with("aggregate_fun_")
  ) %>%
  # omit unnecessary variables
  dplyr::select(
    -data, 
    -province_grid_centers, 
    - fn, 
    -af
  ) %>%
  # change position of the values
  dplyr::relocate(
    ggmap_image_satellite, 
    .after = last_col()
  )
#
# ## save all the results
# Note
# Comment out when not in use. The procedure below takes long computation period.
# # This process needs extremely long periods.
saveRDS(th_sf_rgb, "th_sf_rgb.rds")
gc()
gc()
gc()
#

# Rasterise the obtained satellite imagery

th_sf_rgb <- readRDS("th_sf_rgb.rds")
#
th_sf_raster <-
  th_sf_rgb %>%
  dplyr::mutate(
    image_raster = purrr::map(
      ggmap_image_satellite,
      function(map){
        map_bbox <-
          # Add "boundary box" to the ggmap objects
          base::attr(
            map,
            "bb"
          )
        # change range of the boundary box of ggmap object
        .extent <-
          raster::extent(
            as.numeric(map_bbox[c(2,4,1,3)])
          )
        # rasterise the ggmap object
        my_map <-
          raster::raster(
            x = .extent,
            nrow= nrow(map),
            ncol = ncol(map)
          )
        # obtain color information by RGB format
        rgb_cols <- stats::setNames(as.data.frame(t(grDevices::col2rgb(map))), c('red','green','blue'))
        red <- my_map
        raster::values(red) <- rgb_cols[["red"]]
        green <- my_map
        values(green) <- rgb_cols[["green"]]
        blue <- my_map
        values(blue) <- rgb_cols[["blue"]]
        stack(red,green,blue)
      }
    )
  ) %>%
  # aggregate the imagery using raster:;aggregate() function
  dplyr::mutate(
    aggregate_image_raster = purrr::pmap(
      list(
        image_raster = image_raster,
        aggregate_fun=aggregate_fun,
        fact_number=fact_number
      ),
      function(image_raster, aggregate_fun, fact_number){
        raster::aggregate(x = image_raster,
                          fun = aggregate_fun,
                          fact = fact_number
        )
      }
    )
  )
# ## save all the results
# Note
# Comment out when not in use. The procedure below takes long computation period.
# # This process needs extremely long periods.
saveRDS(th_sf_raster, "th_sf_raster.rds")
gc()
gc()
gc()
#

## make a dataset

th_sf_raster <- readRDS("th_sf_raster.rds")
gc()
gc()
gc()
#
# # Magic spelling for multisession computation
# future::plan(multisession, workers = 4)
# plan()
#
# Transform the rasterized imagery into point with RGB numbers data
th_sf_data <-
  th_sf_raster %>%
  dplyr::select(
    NAME_1,
    individual,
    fact_number,
    aggregate_fun,
    aggregate_image_raster
  ) %>%
  mutate(
    aggregate_image_rasterpoint = furrr::future_map(
      aggregate_image_raster,
      ~
        as_tibble(
          rasterToPoints(.)
        )
    )
  ) %>%
  # omit unnecessary valuables, nest, and groups
  tidyr::unnest(aggregate_image_rasterpoint) %>%
  dplyr::select(-aggregate_image_raster) %>%
  ungroup(.)
#
# ## save all the results
# Note
# Comment out when not in use. The procedure below takes long computation period.
# This process needs extremely long periods.
saveRDS(th_sf_data, "th_sf_data.rds")
# garbage collection
gc()
gc()
gc()
#

# ## make rasterplots from the data above for confirmation

# read original data
th_sf_data <- readRDS("th_sf_data.rds") 
# compute the index
# NOTE
# So many kinds of index are there!!
th_sf_image_index <- 
  th_sf_data %>%
  # limit cases just for trial
  # Comment out when not in use.
  # dplyr::filter(fact_number == 40) %>% 
  # replace variables' name
  data.table::setnames(
    c(
      "province",
      "individual",
      "fact_number",
      "aggregate_fun",
      "lon",
      "lat",
      "red",
      "green",
      "blue"
      )
    ) %>% 
  # # to avoid dividion by zero, add 1 to teh RGB element
  # dplyr::mutate(
  #   red = red + 1,
  #   green = green + 1,
  #   blue = blue + 1,
  # ) %>% 
  group_by(province, individual) %>% 
  nest() %>% 
  # limit cases just for trial
  # Comment out when not in use.
  # dplyr::filter(individual == "1") %>% 
  dplyr::mutate(
    scale.red = purrr::map(
      data,
      ~
        scale(.$red)
    ),
    scale.green = purrr::map(
      data,
      ~
        scale(.$green)
    ),  
    scale.blue = purrr::map(
      data,
      ~
        scale(.$blue)
    )
  ) %>% 
  ungroup() %>% 
  unnest(
    cols = c(
      "data",
      "scale.red",
      "scale.green",
      "scale.blue"
      )
    ) %>% 
  data.table::setnames(
    c(
      "province",
      "individual",
      "fact_number",
      "aggregate_fun",
      "lon",
      "lat",
      "red",
      "green",
      "blue",
      "scale.red",
      "scale.green",
      "scale.blue"
    )
  ) %>% 
  dplyr::mutate(
    red = round(as.numeric(red)),
    green = round(as.numeric(green)),
    blue = round(as.numeric(blue)),
    scale.red = as.numeric(scale.red),
    scale.green = as.numeric(scale.green),
    scale.blue = as.numeric(scale.blue)
  ) %>% 
  # Warning!!
  # To avoid numerical calculation error, we should round the numerator / denominator both
  # In detail, refer to the following page.
  # http://cse.naro.affrc.go.jp/takezawa/r-tips/r/37.html
  dplyr::mutate(
    vari = round((green - red)) / round((green + red - blue)),
    scale.vari = round((scale.green - scale.red) / (scale.green + scale.red - scale.blue)) 
  ) %>% 
  # No matter how we make effort to avoid divergence, it often emerges.
  # To omit observation including that, we should detect that and omit. 
  dplyr::mutate(
    fugainf = is.finite(vari),
    fuganan = is.nan(vari)
  ) %>% 
  dplyr::filter(fugainf == "TRUE" & fuganan == "FALSE") %>% 
  dplyr::select(-fugainf, -fuganan)
# # ## save all the results
# # Note
# # Comment out when not in use. The procedure below takes long computation period.
# # # This process needs extremely long periods.
# saveRDS(th_sf_image_index, "th_sf_image_index.rds")
# gc()
# gc()
# gc()
# 
# #
# ##
# ### --- END --- ###