---
title: "Untitled"
format: html
editor: visual
---

```{r}
#install.packages("stars")
#install.packages("tmap")
#remotes::install_github("r-tmap/tmaptools")
#remotes::install_github("r-tmap/tmap")

```

```{r}
library(rstac)
library(gdalcubes)
library(stars)
library(tmap)
library(dplyr)
```

```{r}
url <- "/vsicurl/https://dsl.richmond.edu/panorama/redlining/static/mappinginequality.gpkg" #places (polygons) from redlining done in the great depression. #vsicurl = virtual file system interface; will read over the command line url and pretend that the 
redlines <- sf::st_read(url) |> sf::st_make_valid() |> filter(st_is_valid(geom))
sf_redlines <- redlines |> filter(city == "San Francisco")
la_redlines <- redlines |> filter(city == "Los Angeles")
```

```{r}
box <- st_bbox(sf_redlines)
start_date <- "2022-06-01"
end_date <- "2022-07-01"
items <-
  stac("https://earth-search.aws.element84.com/v0/") |>
  stac_search(collections = "sentinel-s2-l2a-cogs",
              bbox = c(box),
              datetime = paste(start_date, end_date, sep="/"),
              limit = 100) |>
  ext_query("eo:cloud_cover" < 20) |>
  post_request()
```

```{r}
col <- stac_image_collection(items$features, asset_names = c("B08", "B04", "SCL"))

cube <- cube_view(srs ="EPSG:4326", #Desired datacube
                  extent = list(t0 = start_date, t1 = end_date,
                                left = box[1], right = box[3],
                                top = box[4], bottom = box[2]),
                  dx = 0.0001, dy = 0.0001, dt = "P1D",
                  aggregation = "median", resampling = "average")

mask <- image_mask("SCL", values=c(3, 8, 9)) # mask clouds and cloud shadows

data <-  raster_cube(col, cube, mask = mask)
```

```{r}
ndvi <- data |>
  select_bands(c("B04", "B08")) |>
  apply_pixel("(B08-B04)/(B08+B04)", "NDVI") |>
  reduce_time(c("mean(NDVI)"))

ndvi_stars <- st_as_stars(ndvi)
```

```{r}
mako <- tm_scale_continuous(values = viridisLite::mako(30))
fill <- tm_scale_continuous(values = "Greens")

tmap_mode("plot") #order of the layers matters!
tm_shape(ndvi_stars) + tm_raster(col.scale = mako) + tm_shape(sf_redlines) + tm_borders("grade", lwd=2)
```

```{r}
tmap_mode("view")
tm_basemap() + tm_shape(sf_redlines) + tm_polygons("grade", alpha = 0.5) #tmap() layers maps using ggplot-like syntax
```

```{r}
ok <- redlines |> filter(st_is_valid(geom))
```

```{r}
box <- c(xmin=-122.51, ymin=37.71, xmax=-122.36, ymax=37.81) 
start_date <- "2022-06-01"
end_date <- "2022-08-01"
items <-
  stac("https://earth-search.aws.element84.com/v0/") |>
  stac_search(collections = "sentinel-s2-l2a-cogs",
              bbox = box,
              datetime = paste(start_date, end_date, sep="/"),
              limit = 100) |>
  ext_query("eo:cloud_cover" < 20) |>
  post_request()
```

```{r}
col <- stac_image_collection(items$features, asset_names = c("B08", "B04", "SCL"))

cube <- cube_view(srs ="EPSG:4326",
                  extent = list(t0 = start_date, t1 = end_date,
                                left = box[1], right = box[3],
                                top = box[4], bottom = box[2]),
                  dx = 0.0001, dy = 0.0001, dt = "P1D",
                  aggregation = "median", resampling = "average")

mask <- image_mask("SCL", values=c(3, 8, 9)) # mask clouds and cloud shadows

data <-  raster_cube(col, cube, mask = mask)
```

```{r}
ndvi <- data |>
  select_bands(c("B04", "B08")) |>
  apply_pixel("(B08-B04)/(B08+B04)", "NDVI") |>
  reduce_time(c("mean(NDVI)"))
```
