README.Rmd

# Green Exposure

<!-- badges: start -->
<!-- badges: end -->
- [Installation](#installation)
- [Data](#data)
  * [Ams_Neighborhoods](#ams_neighborhoods)
  * [Ams_Houses](#ams_houses)
  * [Ams_greenspaces](#ams_greenspaces)
- [Functionalities](#functionalities)
  * [Preparation](#preparation)
  * [Availability](#availability)
    + [Calc NDVI](#calc-ndvi)
    + [Land Cover](#land-cover)
    + [Canopy coverage](#canopy-coverage)
    + [greenspace percentage](#greenspace-percentage)
  * [Accessibility](#accessibility)
    + [greenspace access](#greenspace-access)
  * [Visibility](#visibility)
    + [Viewshed](#viewshed)
    + [VGVI](#vgvi)
    + [Streetview](#streetview)
- [Alternative Installation](#alternative-installation)
  * [GEE](#gee)
  * [Rcpp](#rcpp)


# Installation

You can install the development version of GreenExp from [GitHub](https://github.com/) with:

``` r
# install.packages("devtools")
devtools::install_github("Spatial-Data-Science-and-GEO-AI-Lab/GreenEx_R")
```


---

# Data 

This package is provided with three sf datasets:

1. [Ams_Neighborhoods](#ams_neighborhoods)
2. [Ams_Houses](#ams_houses)
3. [Ams_greenspaces](#ams_greenspaces)

These datasets will be used as example to explain the [Functionalities](#functionalities).

#### Ams_Neighborhoods

The first dataset is an sf data frame of neighborhoods in Amsterdam. This dataset is retrieved from [Gemeente Amsterdam](https://maps.amsterdam.nl/gebiedsindeling/). 

Run the following code for more information 

```{r}
library(GreenExp) # If GreenExp is not loaded yet
?Ams_Neighborhoods
```


---

#### Ams_Houses

The Ams_Houses contain houses that are created by taking the centroid of the aforementioned [Ams_Neighborhoods](#Ams_Neighborhoods)
Run the following code for more information 
```{r}
?Ams_Houses
```

---

#### Ams_greenspaces

The Ams_greenspaces dataset is also retrieved from the [Gemeente Amsterdam](https://maps.amsterdam.nl/stadsgreenspaceen/)
Run the following code for more information 

```{r}
?Ams_greenspaces
```

---


# Functionalities 

---

## Preparation

The Functionalities, which will be treated in the next subsections, will be provided with examples.
To avoid computationally heavy examples, a few neighborhoods in Amsterdam will be selected, namely: 
Rapenburg, Uilenburg, Valkenburg, Marine-Etablissement, Kadijken, Plantage, Kattenburg, Wittenburg and Oostenburg. 

The code below makes this dataset and showcases in the plot which neighborhoods are selected.

```{r, warning =FALSE, message=FALSE}
library(GreenExp) # If not loaded yet
library(magrittr) # If not loaded yet (used for piping %>%)


neighborhoods <- c('Rapenburg', 'Uilenburg', 'Valkenburg', 
            'Marine-Etablissement', 'Kadijken', 'Plantage', 
            'Kattenburg', 'Wittenburg', 'Oostenburg')

# Filter the neighborhoods
df <- Ams_Neighborhoods %>%
  dplyr::filter(Buurt %in% neighborhoods)

# Create point locations in Amsterdam  
df_points <- sf::st_centroid(df)

# Create the greenspaces
#df_greenspaces <- Ams_greenspaces



map <- mapview::mapview(list(df, df_points), label = c("Buurt", "Buurtcode"),
                        alpha.regions=c(0,1),  legend=F, lwd=1)

# Display the updated map
map
```

Please note that the data serves as an illustration, and you may need to adapt the parameters and function usage to match your specific scenario.

---

## Availability

Availability will be assessed using four functions:

1. [Calc NDVI](#calc-ndvi)
2. [Land Cover](#land-cover)
3. [Canopy coverage](#canopy-coverage)
4. [greenspace percentage](#greenspace-percentage) 

Each of these functions will provide an [sf](https://r-spatial.github.io/sf/articles/sf1.html) `dataframe` that includes the input location and the specific values requested within a defined buffer.


The user has the option to input either a point geometry or a (multi)polygon geometry. By default, the address location will be transformed into a point geometry, and a buffer will be created around it to calculate the availability. However, users can choose to use the provided polygon geometry to calculate availability by setting the 'address_location_neighborhood' parameter to TRUE.

By default, the buffer around the input location is measured in Euclidean distance. However, it can be modified to utilize a network buffer. The distinction between the two types of buffers is illustrated in the figure below. The Euclidean buffer in this instance has a fixed radius of 1000 meters, while the network buffer is calculated based on a speed of 5 km/h over a duration of 10 minutes.

![](man/figures/buffers_example.png)


In the following subsections, a brief description of each availability function will be provided, along with examples extracted from the neighborhood polygons and points in Amsterdam. 

---

### Calc NDVI 

The `calc_ndvi` function computes the average Normalized Difference Vegetation Index [(NDVI)](https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index) within a specified distance for given location(s). The input for the function is `address_location` which should be an `sf dataframe`. It is recommended to provide the `address location` with a projected Coordinate Reference System [(CRS)](https://docs.qgis.org/3.28/en/docs/gentle_gis_introduction/coordinate_reference_systems.html#:~:text=In%20layman%27s%20term%2C%20map%20projections,real%20places%20on%20the%20earth). If no projected CRS is provided, the address location will be automatically projected to [WGS 84 / World Mercator](https://epsg.io/3395). 


You have the option to provide a raster file containing NDVI values. However, if no raster file is provided, the function will use the [Sentinel-2-l2a](https://planetarycomputer.microsoft.com/dataset/sentinel-2-l2a) dataset from Planetary Computer as the default data source for calculating NDVI. The figures below illustrate two examples of NDVI in Amsterdam. The first plot showcases the aforementioned neighborhoods in Amsterdam and there NDVI, the second one shows the address points within an euclidean buffer of 300 meters 

Neighborhood NDVI        |  Addresses within euclidean distance NDVI
:-------------------------:|:-------------------------:
![](man/figures/NDVI_neighborhoods.png)   |  ![](man/figures/NDVI_points_euclidean.png)

If desired, users are able to switch the engine to `GEE` (Google Earth Engine) for performing the calculations


---

Below you can find the code where the results correspond with the NDVI figure. 

```{r, warning =FALSE, message=FALSE}
ndvi_neighborhood <- GreenExp::calc_ndvi(df, address_location_neighborhood = TRUE, UID=df$Buurt)
ndvi_neighborhood
```

```{r, warning =FALSE, message=FALSE}
ndvi_address <- GreenExp::calc_ndvi(df_points, buffer_distance=300)
ndvi_address
```

---

### Land Cover


The `land_cover` function calculates the average land cover within a specified distance for given location(s). The input for the function is address_location, which should be an `sf dataframe` with projected Coordinate Reference System (CRS) information. If no CRS is provided, the function will automatically project the address location to WGS 84 / World Mercator.

You have the option to provide a raster file wiht land cover values. When this is not provided, the [esa-worldcover](https://planetarycomputer.microsoft.com/dataset/esa-worldcover) data set of Planetary Computer will be used to calculate the land cover.

In the code chunk and figure below an example is given for the Amsterdam area. It illustrates an example of land cover within a network buffer of 500 meters in Amsterdam. It showcases three address locations, and the land cover is determined based on the surrounding road network obtained from the  [osmextract](https://cran.r-project.org/web/packages/osmextract/vignettes/osmextract.html) package.


![](man/figures/land_cover_network.png) 

---
```{r, warning=FALSE, message=FALSE}
network_file <- network_file
land_cover_values <- GreenExp::land_cover(df_points, buffer_distance=300, network_buffer=T, network_file = network_file)
land_cover_values 
```


---

### Canopy coverage

The `canopy_perc` function calculates the percentage of a canopy within a given buffer distance or location. 

``` r
# Read the canopy dataset
canopy <- sf::st_read("Data/CanopyTestArea.gpkg")

canopy_perc(address_location = address_test, canopy_layer = canopy, buffer_distance = 500)
 

# Simple feature collection with 3 features and 2 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 385981.9 ymin: 392861.6 xmax: 388644.2 ymax: 395322.2
# Projected CRS: OSGB36 / British National Grid
#   UID canopy_pct                  geometry
# 1   1   14.42063 POINT (388644.2 392861.6)
# 2   2   19.27852 POINT (385981.9 393805.5)
# 3   3   10.67145 POINT (388631.2 395322.2)
```

---

### greenspace percentage

The `greenspace_pct` function gives the percentage of greenspace coverage given a certain buffer. If the `greenspace_layer` is not given, the greenspaces will be retrieved using features from [osmdata](https://wiki.openstreetmap.org/wiki/Map_features). 

![](man/figures/parks_euclidean.png)


```{r, warning =FALSE, message=FALSE}
greenspace_pct(df_points, speed=5, time=3)
```


---

## Accessibility
### Greenspace access


The greenspace_access function provide the ability to determine the nearest greenspaces to given address locations and assess their accessibility within a specified buffer distance. By default, the functions utilize a euclidean buffer around the address locations and calculate the shortest distance to the centroid of the greenspaces. This is achieved using the K-nearest neighbors (KNN) algorithm with the [FNN](https://rdrr.io/cran/FNN/man/knn.html) package, to calculate the euclidean distance between the address location and the greenspaces.


Furthermore, the functions allow for the option to utilize a network buffer instead of the euclidean buffer. In this case, the distance calculation is performed using the [sfnetworks](https://cran.r-project.org/web/packages/sfnetworks/index.html) package, which leverages road networks to calculate distances between points.

Additionally, pseudo entry points can be employed to calculate the distance to the greenspaces. These pseudo entrances are created by generating a 10-meter buffer around the greenspace polygons and intersecting them with the network nodes obtained from the intersection of the network points with the greenspaces.


Three examples will be provided. For the sake of visualization, one address point will be used to calculate the accessibility. 

**Example 1: Euclidean Distance Calculation**

In this example, the accessibility function is applied using the default settings, which involves calculating the euclidean distance from the address location to the nearest greenspace centroid. The figure below illustrates an example in Amsterdam, where the greenspaces are represented by green polygons. The blue lines indicate the euclidean distance from the address location to the nearest greenspace centroid. The greenspace centroids are depicted as black points, while the address location is denoted by a red point. The code chunk beneath the plot provides the necessary code to receive the shortest distance from the address location. 



```{r, include=FALSE}
library(sf)
library(mapview)
df_point <- df_points[8,]

projected_crs <- sf::st_crs(df_point)

address_locations <- df_point

buffer_distance <- 300

calculation_area <- sf::st_buffer(address_locations, dist = buffer_distance)

calculation_area <- sf::st_transform(calculation_area, 4326)

iso_area <- sf::st_buffer(sf::st_convex_hull(
    sf::st_union(sf::st_geometry(df_point))),
    buffer_distance)
  iso_area <- sf::st_transform(iso_area, crs = 4326)
q1 <- osmdata::opq(sf::st_bbox(iso_area)) %>%
      osmdata::add_osm_feature(key = "landuse",
                               value = c('allotments','forest',
                                         'greenfield','village_green')) %>%
      osmdata::osmdata_sf()

q2 <- osmdata::opq(sf::st_bbox(iso_area)) %>%
  osmdata::add_osm_feature(key = "leisure",
                           value = c('garden','fitness_station',
                                     'nature_reserve','park','playground')) %>%
  osmdata::osmdata_sf()

q3 <- osmdata::opq(sf::st_bbox(iso_area)) %>%
  osmdata::add_osm_feature(key = "natural",
                           value = c('grassland')) %>%
  osmdata::osmdata_sf()
res <- c(q1, q2, q3)



greenspace <- res$osm_polygons
greenspace <- tidygraph::select(greenspace, "osm_id", "name")
greenspace <- sf::st_make_valid(greenspace)

greenspace_centroid <- sf::st_centroid(greenspace)
address <- sf::st_transform(address_locations, 4326)

# Calculate Euclidean distances between address and greenspace centroids
# Multiple distances 
nearest_neighbours <- FNN::get.knnx(sf::st_coordinates(greenspace_centroid), sf::st_coordinates(address), k = 5)

indices <- nearest_neighbours$nn.index[, 1]  # Extract the indices from the first c

greenspace_centroids <- sf::st_as_sf(greenspace_centroid)

greenspace_centroid_matrix <- as.matrix(greenspace_centroid)
selected_rows <- greenspace_centroid_matrix[indices, , drop = FALSE]



distances <- st_distance(address, greenspace_centroid)

# Find the index of the greenspace centroid with minimum distance
closest_index <- which.min(distances)

greenspace_centroid <- sf::st_as_sf(greenspace_centroid)
closest_greenspace <- greenspace_centroid[closest_index,]

calc_area <- sf::st_transform(calculation_area, projected_crs)
greenspace_centroid <- sf::st_transform(greenspace_centroid, projected_crs)
closest_greenspace <- sf::st_transform(closest_greenspace, projected_crs)

coords_address <- data.frame(sf::st_coordinates(address_locations))
closest_greenspace_coords <- data.frame(sf::st_coordinates(closest_greenspace))

line <- st_linestring(matrix(c(coords_address$X, closest_greenspace_coords$X, coords_address$Y, closest_greenspace_coords$Y), ncol = 2))
line <- st_sfc(line, crs = projected_crs)
line <- sf::st_as_sf(line)

greenspace_layer <- sf::st_transform(greenspace, projected_crs)
# Display the mapview plot with the legend
m1 <- mapview(greenspace_layer, col.regions='green', layer.name='Greenspace Polygon', legend=TRUE)
m2 <- mapview(greenspace_centroid, col.regions='black', cex=3, layer.name='Greenspace Centroid')
m3 <- mapview(address_locations, col.regions='red', cex=3, layer.name='Address Location')
m4 <- mapview(line, layer.name='Euclidean line')

map <- m1+m2+m3+m4

```

```{r}
map
```

```{r, warning =FALSE, message=FALSE}
greenspace_access(df_point, buffer_distance = 300)
```

---

**Example 2: Network Distance Calculation**

In this example, the accessibility function utilizes network distance to compute the distance from the address location to the nearest greenspace centroid. The figure below showcases an example in Amsterdam, where the greenspaces are represented by green polygons. However, since the lines are retrieved from an OSM network file, the greenspace centroids and address centroid may not align exactly with the network lines. As a result, you may notice that the lines do not intersect with the points in the plot. The greenspace centroids are depicted as black points. The code chunk below the plot corresponds with the distance result of the plot




```{r, include=FALSE}
lines <- GreenExp::network_file
iso_area <- sf::st_buffer(sf::st_convex_hull(
    sf::st_union(sf::st_geometry(df_point))),
    buffer_distance)
  iso_area <- sf::st_transform(iso_area, crs = 4326)
region_shp <- sf::st_transform(iso_area, projected_crs)

#Download osm used a square bounding box, now trim to the exact boundary
#note that lines that that cross the boundary are still included
lines <- lines[region_shp,]

# Round coordinates to 0 digits.
sf::st_geometry(lines) <- sf::st_geometry(lines) %>%
  sf::st_sfc(crs = sf::st_crs(lines))

# Network
network_file <- sfnetworks::as_sfnetwork(lines, directed = FALSE)
network_file <- tidygraph::convert(network_file, sfnetworks::to_spatial_subdivision)

#convert network to an sf object of edge

net_sf <- network_file %>% tidygraph::activate("edges") %>%
  sf::st_as_sf()

# Find which edges are touching each other
touching_list <- sf::st_touches(net_sf)
# create a graph from the touching list
graph_list <- igraph::graph.adjlist(touching_list)
# Identify the cpnnected components of the graph
roads_group <- igraph::components(graph_list)
# cont the number of edges in each component
roads_table <- table(roads_group$membership)
#order the components by size, largest to smallest
roads_table_order <- roads_table[order(roads_table, decreasing = TRUE)]
# get the name of the largest component
biggest_group <- names(roads_table_order[1])

# Subset the edges corresponding to the biggest connected component
osm_connected_edges <- net_sf[roads_group$membership == biggest_group, ]
# Filter nodes that are not connected to the biggest connected component
network_file <- network_file %>%
  tidygraph::activate("nodes") %>%
  sf::st_filter(osm_connected_edges, .pred = sf::st_intersects)

net_intersect <- sf::st_transform(network_file,  crs = 4326)
net_points <- sf::st_geometry(net_intersect)
# Compute the edge weights bsased on their length
network_file <- tidygraph::mutate(tidygraph::activate(network_file, "edges"),
                                  weight = sfnetworks::edge_length())

distance_greenspaces<-sfnetworks::st_network_cost(network_file, from = df_point, to = greenspace_centroid)
shortest_distances <- apply(distance_greenspaces, 1, which.min)
greenspace_centroid_closest <- greenspace_centroid[shortest_distances, ]


closest_greenspace <- sfnetworks::st_network_paths(network_file, from=df_point, to=greenspace_centroid_closest)
node_path <- closest_greenspace  %>%
  dplyr::slice(1) %>%
  dplyr::pull(node_paths) %>%
  unlist()
line_network <- dplyr::slice(tidygraph::activate(network_file, "nodes"), node_path)
final_line<- line_network %>% tidygraph::activate("edges") %>%
       sf::st_as_sf() 
m4 <- mapview(final_line, layer.name='Network line')

map <- m1+m2+m3+m4
```

```{r}
map
```

```{r, warning =FALSE, message=FALSE}
greenspace_access(df_point, buffer_distance = 300, euclidean=F, network_file = GreenExp::network_file)
```

---

**Example 3: Network Distance to Pseudo Entrances**

In this example, the accessibility function considers network distance to the pseudo entrances of the greenspaces. The pseudo entrances are created by buffering the greenspace polygons and intersecting them with the network nodes. The function calculates the network distance from the address location to the nearest pseudo entrance point. The figure below presents an example in Amsterdam, where the greenspaces are shown as green polygons. The blue lines indicate the euclidean distance from the address location to the nearest greenspace centroid. The greenspace centroids are depicted as black points, and the address location is represented by a red point. Additionally, you may observe multiple pseudo entrances within the greenspaces, as roads passing through the greenspaces can also serve as potential entrance points. 

```{r, include=FALSE}
greenspace_buffer <- sf::st_buffer(greenspace, 20)

# When greenspace are overlapping, combine them
greenspace_combined <- sf::st_union(greenspace_buffer)

# make it a multipolygon
greenspace_buffer <- sf::st_as_sf(sf::st_sfc(greenspace_combined))
greenspace_buffer <- sf::st_make_valid(greenspace_buffer)

st_intersection_faster <- function(x,y){
  #faster replacement for st_intersection(x, y,...)
  
  y_subset <-
    sf::st_intersects(x, y) %>%
    unlist() %>%
    unique() %>%
    sort() %>%
    {y[.,]}
  
  sf::st_intersection(x, y_subset)
}
greenspace_point <- st_intersection_faster(net_points, greenspace_buffer)
greenspace_point <- greenspace_point %>% sf::st_transform(crs=projected_crs)


distance_greenspaces<-sfnetworks::st_network_cost(network_file, from = df_point, to = greenspace_point)
shortest_distances <- apply(distance_greenspaces, 1, which.min)
greenspace_point <- sf::st_as_sf(greenspace_point)
greenspace_centroid_closest <- greenspace_point[shortest_distances, ]


closest_greenspace <- sfnetworks::st_network_paths(network_file, from=df_point, to=greenspace_centroid_closest)
node_path <- closest_greenspace  %>%
  dplyr::slice(1) %>%
  dplyr::pull(node_paths) %>%
  unlist()
line_network <- dplyr::slice(tidygraph::activate(network_file, "nodes"), node_path)
final_line<- line_network %>% tidygraph::activate("edges") %>%
       sf::st_as_sf() 
m4 <- mapview(final_line, layer.name='Network line')


greenspace_pseudo <- sf::st_as_sf(greenspace_point)

m2 <- mapview(greenspace_pseudo, col.regions='black', cex=3, layer.name='Pseudo entrances')

map <- m1 + m2 + m3 + m4
```

```{r}
map
```



```{r, warning =FALSE, message=FALSE}
greenspace_access(df_point, buffer_distance=300, euclidean = F, pseudo_entrance = T, 
                  network_file = GreenExp::network_file)
```


## Visibility
 
The visibility functions are made by the [GVI](https://github.com/STBrinkmann/GVI) package with some adaptations. 

---

### Viewshed

The viewshed function computes a binary viewshed of a single point on a Digital Surface Model (DSM) raster. A radial buffer is applied on the observer position, and visibility is being calculated using a C++ implementation of Bresenham’s line algorithm [Bresenham 1965](https://ieeexplore.ieee.org/document/5388473) & [Bresenham 1977](https://doi.org/10.1145/359423.359432) and simple geometry. The
result of the `viewshed` function is a radial raster where 0 = no-visible and 1 = visible area.

For a better explanation, go to the [GVI](https://github.com/STBrinkmann/GVI) package.

**EXAMPLE**

```{r}
# Read the greenspace 
GS <- terra::rast('data/GS_AMS.tif')
# Read the digital surfaca model
DSM <- terra::rast('data/DSM_AMS.tif')
# Read in the digital eleveation model
DTM<- terra::rast('data/DTM_AMS.tif')
vs <- GreenExp::viewshed(observer = df_points[1,], dsm_rast = DSM, dtm_rast = DTM,
                         max_distance = 200, observer_height = 1.7, plot = TRUE)
```


The left plot represents the Digital Surface Model (DSM), whereas the right plot represents the viewshed, where green is the visibile area and gray is not visible. 

### VGVI

The Viewshed Greenness Visibility Index (VGVI) represents the proportion of visible greenness to the total visible area based on the `viewshed`. The estimated VGVI values range between 0 and 1, where = no green cells and 1= all of the visible cells are green.

Based on a viewshed and a binary greenspace raster, all visible points are classified as visible green and visible no-green. All values are summarized using a decay function, to account for the reducing visual prominence of an object in space with increasing distance from the observer. Currently two options are supported, a logistic and an exponential function.

For more information about the VGVI please go to the [GVI](https://github.com/STBrinkmann/GVI) package. For more information about the algorithms look at the paper by [Brinkmann, 2022](https://doi.org/10.5194/agile-giss-3-27-2022)

**Example** 

The figure below provides an overview of the data used for calculating the Vegetation-Ground View Index (VGVI). In the figure, you will find various elements that contribute to the analysis.

1. Observer Points: The dots depicted in the figure represent the observer points. These points correspond to the addresses used in the examples discussed thus far. Each observer point serves as a starting location for measuring the VGVI.

2. Address ID Numbers: The numbers assigned to the addresses in the plot correspond to the ID numbers used in the code chunk below the figure. These ID numbers uniquely identify each address and allow for easy referencing and analysis in the code.

3. Green Space: The green shades in the plot represent the tree raster, which represents the extent of green space. It indicates the areas covered by vegetation, such as trees or greenspaces, and is a crucial factor in determining the VGVI.DTM and 

4. DSM: The figure also includes the Digital Terrain Model (DTM) and Digital Surface Model (DSM). These models provide information about the elevation of the terrain and structures present in the area. The combination of DTM and DSM helps in understanding the topography of the region.

By utilizing this information and the corresponding code, the VGVI can be calculated, providing insights into the vegetation-ground view characteristics at each observer point.




![](man/figures/VGVI_AMS.png)

```{r}
VGVI <- GreenExp::vgvi_from_sf(observer = df_points,
                     dsm_rast = DSM, dtm_rast = DTM, greenspace_rast = GS,
                     max_distance = 200, observer_height = 1.7,
                     m = 0.5, b = 8, mode = "logit")

VGVI

```



--- 

### VGVI from address

GVI from address function: In contrast, the VGVI from address function employs a broader approach. It samples multiple points around the address location within a defined buffer. This buffer represents a circular area around the address. The function collects data from various points within this buffer and calculates the VGVI by taking the mean of the collected values. By incorporating multiple sample points, it offers a more comprehensive representation of the VGVI within the vicinity of the address.

The VGVI from sf function analyzes the VGVI at a specific observer point, while the VGVI from address function expands the analysis by sampling multiple points around the address location. The latter approach captures the variability in vegetation-ground view characteristics within a defined buffer and provides a more averaged assessment of the VGVI.


```{r}
mean_VGVI <- vgvi_from_address(address = df_points,
                     dsm_rast = DSM, dtm_rast = DTM, greenspace_rast = GS,
                     max_distance = 200, observer_height = 1.7,
                     m = 0.5, b = 8, mode = "logit")

mean_VGVI
```

---

# Alternative installation

To make optimal use of the package 

## GEE 

This step is optional, by default the [Planetary Computer](https://planetarycomputer.microsoft.com) will be used for satellite images, 
But if you also want to use the [Google Earth Engine](https://earthengine.google.com), and do not have it installed yet,
you need to follow the following steps or the steps given in this [instruction video](https://www.youtube.com/watch?v=_fDhRL_LBdQ)

**Step 1:**

Make an account on  [Google Earth Engine](https://earthengine.google.com)


``` r
install.packages(c("sf", "reticulate", "remotes"))
```
afterwards install the [rgee](https://github.com/r-spatial/rgee) package from github

``` r
# Install the rgee package from the r spatial GitHub
remotes::install_github("r-spatial/rgee")

# load the reticulate and rgee package
library(reticulate)
library(rgee)
```

**Step 2:**

Running `reticulate::py_discover_config()` will install `Miniconda`

``` r
# Use the py_discover_config() function to see what version of Python will be used
# without actually loading pythong
reticulate::py_discover_config()

# python:         /Users/martijn/.virtualenvs/rgee/bin/python
# libpython:      /Users/martijn/anaconda3/lib/libpython3.10.dylib
# pythonhome:     /Users/martijn/.virtualenvs/rgee:/Users/martijn/.virtualenvs/rgee
# version:        3.10.9 (main, Mar  1 2023, 12:20:14) [Clang 14.0.6 ]
# numpy:          /Users/martijn/.virtualenvs/rgee/lib/python3.10/site-packages/numpy
# numpy_version:  1.24.3

# Verify the current Python path
import('sys')$executable
# [1] "/Users/martijn/.virtualenvs/rgee/bin/python"


# Create an isolated Python venv with all rgee dependencies
ee_install()
# look at the path to the rgee env
```

**Step 3:**

After this bit, please restart your pc/laptop and launch R again. 

**Step 4:**

Initializing

``` r
# Set python version to use
reticulate::use_python("/Users/martijn/.virtualenvs/rgee/bin/python")
reticulate::py_config()

library(rgee)

#Initialize the Earth Engine
ee_Initialize()

## 2. Install geemap in the same Python ENV that use rgee
py_install("geemap")
gm <- import("geemap")

```

Enter the email you used to sign-up for GEE 

copy the code into R 

── rgee 1.1.6.9999 ───────────────────────────────────── earthengine-api 0.1.354 ── 

 ✔ user: not_defined 
 
 ✔ Initializing Google Earth Engine:  DONE!
 
 ✔ Earth Engine account: users/ee-greenexp 
 
---

## Rcpp 

Make sure the [Rcpp package](https://cran.r-project.org/web/packages/Rcpp/index.html) is installed.
If you are using mac, make sure you have [Xcode](https://apps.apple.com/nl/app/xcode/id497799835?mt=12) installed. 

Furthermore you have to make a Makevars file if the cpp files are not working. 
go to terminal and do the following:

```
mkdir .R
cd .R 
touch Makevars
open Makevars

## copy and paste:
FC = /opt/homebrew/Cellar/gcc/13.1.0/bin/gfortran
F77 = /opt/homebrew/Cellar/gcc/13.1.0/bin/gfortran
FLIBS = -L/opt/homebrew/Cellar/gcc/13.1.0/lib/gcc/13
```
---

## Mapillary

To use the [streetview](#streetview) function, data will be retrieved using the [mapillary](https://www.mapillary.com) API