---
title: "Introduction to greenR"
author: "Sachit Mahajan"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to greenR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# 1. Introduction

`greenR` is an open-source R package designed to quantify, analyze, and visualize urban greenness using data from OpenStreetMap (OSM). This package is essential for urban planners, environmental researchers, and policy-makers interested in urban sustainability and green infrastructure.

This vignette will guide you through the installation and basic usage of `greenR`, providing examples of its main features.

# 2. Installation

You can install the development version of `greenR` from GitHub:

```{r, eval=FALSE}
# install.packages("devtools")
devtools::install_github("sachit27/greenR", dependencies = TRUE)
```

# 3. Download data

The first step is to acquire data. This provides a systematic approach to collecting the requisite geospatial data from OSM, thereby serving as the foundation for all subsequent analyses. The users can simply specify any city or neighborhood (that has data available in OSM database). This function looks in the database and finds any city and downloads OSM data for the specified spatial area with regard to three key environmental features: highways, green areas, and trees. Here green areas include all the areas with the following tags: "forest", "vineyard", "plant_nursery", "orchard", "greenfield", "recreation_ground", "allotments", "meadow","village_green","flowerbed", "grass", "farmland", "garden", "dog_park","nature_reserve", and "park".

```{r, eval = FALSE}
library(greenR)
data <- get_osm_data("City of London, United Kingdom")
```

# 4. Visualize green spaces and clustering

The 'visualize_green_spaces()' function is designed to aid in the visual assessment of green space data. Utilizing an integrated leaflet map, users can explore the distribution and mapping quality of green spaces within a specified area. By plotting this data on an interactive leaflet map, users gain insights into the extent and accuracy of green space representation. After visualizing the green spaces within the desired area, users may wish to contribute to the OpenStreetMap project to enhance the data quality or add unrepresented areas. Additionally, the green_space_clustering() function initially transforms the green spaces into an equal-area projection to calculate the areas accurately. Post-transformation, the K-means algorithm is applied to these areas, clustering the green spaces based on the number of clusters specified.

```{r, eval = FALSE}
green_areas_data <- data$green_areas
visualize_green_spaces(green_areas_data)
green_space_clustering(green_areas_data, num_clusters = 3)
```

# 5. Accessibility analysis

The accessibility_greenspace function creates an interactive map displaying accessible green spaces within a specified walking time from a provided location. It utilizes isochrones to visualize the areas reachable by walking. The function relies on pedestrian routing information and green space data to accurately delineate accessible areas. There are two functions to do this. One with Mapbox and the other with leaflet. For Mapbox, first remember to add the token that you can get from the website.

Mapbox Version (Dynamic): The accessibility_mapbox function creates an accessibility map using Mapbox GL JS, showing green areas and allowing users to generate isochrones for walking times. This produces an HTML that has the inbuilt feature of changing the walking time and also moving the location marker.

```{r, eval = FALSE}
mapbox_token <- "your_mapbox_access_token_here"
accessibility_mapbox(green_areas_data, mapbox_token)
```

Leaflet version: This on euses Leaflet to create the map. In this case, the latitude and longitude has to be given by the user. The default maximum walking time is set to 15 minutes but can be adjusted using the 'max_walk_time' parameter.\ The 'nearest_greenspace' function calculates and visualizes the shortest walking route to the nearest green space of a specified type from a given location. You can also select the type of green space (c("park", "forest")). By default it looks at all the green space.

```{r, eval = FALSE}
accessibility_greenspace(green_areas_data, 47.56427527336772, 7.595820936462059)
map <- nearest_greenspace(data$highways, data$green_areas, 51.4761, -0.2008)
```

```{r, eval = FALSE}
print(map)  # Displays the interactive map
```

# 6. Green Space Accessibility Analysis Using Street Networks

greenR provides a robust framework to assess how easily residents can access green spaces using actual travel routes. This goes beyond isochrones and uses real street networks (for walking, cycling, or driving), as well as optionally population data, to produce scientifically valid, interpretable accessibility metrics.
The recommended workflow is:

-Analyze accessibility using analyze_green_accessibility()

-Create static, radar, and interactive visualizations with create_accessibility_visualizations()

```{r, eval = FALSE}
# Example: Analyze accessibility in the City of London using walking network
library(greenR)
library(terra) # For raster population, if used

# Download OSM data
data <- get_osm_data("City of London, United Kingdom")

# Run accessibility analysis for walking
result <- analyze_green_accessibility(
  network_data = data$highways$osm_lines,
  green_areas = data$green_areas$osm_polygons,
  mode = "walking",
  grid_size = 300 # meters (change as needed)
  # population_raster = my_pop_raster # Optional: a terra::SpatRaster layer
)

# The result is a list with grid-level accessibility and summary stats.
str(result)
```

```{r, eval = FALSE}
# Generate visualizations (static plots, radar plot, interactive Leaflet)
viz <- create_accessibility_visualizations(
  accessibility_analysis = result,
  green_areas = data$green_areas$osm_polygons,
  mode = "walking"
)

# Publication-ready distance map
print(viz$distance_map)

# Barplot: spatial and population-weighted coverage within thresholds
print(viz$coverage_plot)

# Radar plot: directional accessibility (labels inside the axis)
print(viz$directional_plot)

# Combined plot (distance, coverage, radar)
print(viz$combined_plot)

# Interactive map (can be viewed in RStudio Viewer)
viz$leaflet_map

# Summary stats as text
cat(viz$summary)

# Directional table for further analysis
print(viz$directional_table)
```
What you get:

- Static maps of accessibility and coverage
- Directional/radar plot
- Interactive Leaflet map with layers (distance, population, green spaces)
- Textual summary and exportable data

# 7. Calculate the green index for the specified city

This function takes as input the OSM data, a Coordinate Reference System (CRS) code, and parameter D for the distance decay functions. The algorithm extracts the highways, green areas, and trees data from the input list and transforms the data into the given CRS. The CRS affects how distances, areas, and other measurements are calculated. Different CRSs may represent the Earth's surface in ways that either exaggerate or minimize certain dimensions. So, using the wrong CRS can lead to incorrect calculations and analyses. If you're focusing on a city or other localized area, you'll likely want to use a CRS that is tailored to that specific location. This could be a local city grid system or other local CRS that has been designed to minimize distortions in that area. This function then defines distance decay functions for green areas and trees using the parameter D. For each edge in the highway data, the function calculates the green index using the decay functions and returns a data frame with the green index for each edge. By default, D is specified to 100 (distance decay parameter in meters) but it can be changed by the user. Similarly, the users must specify the CRS (<https://epsg.io/>). The green index ranges from 0 to 1 and it represents the relative greenness of each section, factoring in proximity to green spaces and tree density.

```{r, eval = FALSE}

tryCatch({
  green_index <- calculate_green_index(data, 4326, 100)
  print(head(green_index))
}, error = function(e) {
  message("Note: Green index calculation was skipped in this vignette build.")
  message("This function requires the DuckDB spatial extension, which may not be available in all environments.")
  # Create a dummy green_index object for demonstration purposes
  green_index <- data.frame(
    osm_id = 1:5,
    green_index = runif(5),
    geometry = rep("LINESTRING(0 0, 1 1)", 5)
  )
  green_index <- sf::st_as_sf(green_index, wkt = "geometry", crs = 4326)
  print(head(green_index))
})


# 8. Create the green index plot

This function visualizes the green index on a map, with options for both static and interactive display. Interactive maps are rendered using Leaflet, allowing users to zoom, pan, and interact with the map to explore the green index in more detail.

-   Dynamic Mapping: Create interactive, dynamic maps for a more engaging and detailed visualization.
-   Customization: Modify color palette, text size, resolution, title, axis labels, legend position, line width, and line type to suit your preferences.

```{r, eval = FALSE}
# Create a static plot
map <- plot_green_index(green_index)

# Customize static plot
map <- plot_green_index(green_index, colors = c("#FF0000", "#00FF00"), line_width = 1, line_type = "dashed")

# Create an interactive plot using Leaflet
map <- plot_green_index(green_index, interactive = TRUE, base_map = "CartoDB.DarkMatter")

# Use a light-themed base map
map <- plot_green_index(green_index, interactive = TRUE, base_map = "CartoDB.Positron")

```


```{r, eval = FALSE}
print(map) #to view the map in the console. You can use htmlwidgets to save the html
```

# 9. Visualizing data as 3D linestring map

The create_linestring_3D function generates a 3D linestring map using Mapbox GL JS. This map is designed to visualize linear features such as roads, trails, or any other types of linestring data. This can be particularly useful for visualizing connectivity, transportation networks, or other linear spatial patterns.The function supports interactive controls for adjusting the line width and toggling building visibility on the map. It accepts linestring data and automatically processes it to create a visually appealing 3D map.


```{r, eval = FALSE}
mapbox_token <- "your_mapbox_access_token_here"

create_linestring_3D(green_index, "green_index", mapbox_token)
```

# 10. Visualizing green index or any other geospatial data using 3D hexagon map
The 'create_hexmap_3D' function generates a 3D hexagon map using H3 hexagons and Mapbox GL JS. This map can visualize various types of geographical data, such as points, linestrings, polygons, and multipolygons. It is particularly useful for visualizing density or green indices over an area. It automatically processes these geometries, converting them to points for visualization. Users can dynamically change the radius of the hexagons and their heights to better represent the data. The resulting map includes controls for adjusting hexagon height and H3 resolution, and selecting different Mapbox styles.

```{r, eval = FALSE}
mapbox_token <- "your_mapbox_access_token_here"

# Create a temporary file
temp_file <- tempfile(fileext = ".html")

# Create the 3D hexagon map
create_hexmap_3D(
  data = green_index,
  value_col = "green_index",
  mapbox_token = mapbox_token,
  output_file = temp_file,
  color_palette = "interpolateViridis"
)

# Remove the temporary file
unlink(temp_file)
```

# 11. Calculate the percentage of edges with a certain green index

This function groups the edges by their respective green index and calculates the percentage of edges for each green index. For easier interpretation, we categorize the index into three tiers: Low ( < 0.4), Medium (0.4-0.7), and High (> 0.7). 

```{r, eval = FALSE}
percentage <- calculate_percentage(green_index)
```
# 12. Data export and sharing

These functions allow the user to download the green index values as a GeoJson file as well as a Leaflet map. The GeoJSON file retains the geographical properties of the data and can be readily employed in a broad range of GIS applications. The Leaflet map, saved as an HTML file, provides an interactive user experience, facilitating dynamic exploration of the data. The users should specify the file path to save these files.


```{r,eval = FALSE}
download_file <- save_json(green_index, "File_Path") #file path has to be specified. For example "/Users/.../map.geojson"
map <- save_as_leaflet(green_index, "File_Path")
```

# 13. Green and Tree Count Density Analysis (H3 Hexagons, Inequality Metrics, and Maps)

You can assess the spatial distribution and inequality of green spaces or trees using H3 hexagonal binning and get advanced metrics including Lorenz curves and the Gini index. This approach is robust to OSM data gaps and provides publication-ready analytics and visualizations.

```{r, eval = FALSE}
# Example: Analyze green area polygon density
osm_data <- get_osm_data("City of London, United Kingdom")

result <- analyze_green_and_tree_count_density(
  osm_data = osm_data,
  mode = "green_area",      # or "tree_density"
  h3_res = 8,               # H3 resolution (smaller is coarser)
  color_palette = c("#FFEDA0", "#74C476", "#005A32"),
  save_lorenz = TRUE
)

# Interactive Leaflet map (colored by density bins)
result$map

# Print analytics (mean, median, skewness, Gini, count per km2, etc)
print(result$analytics)

# (Optional) Lorenz curve PNG path
result$lorenz_plot
```
You can also analyze tree density by setting mode = "tree_density".

# 14. ALS GEDI Canopy Height Model (CHM) Analysis

greenR enables rapid analysis and visualization of global 1m Canopy Height Model (CHM) data from Meta & WRI's ALS GEDI v6 dataset. Analyze, mosaic, crop, and visualize canopy heights for any city or AOI - including tree cover fraction and AOI size (Mollweide km^2) with maps and histograms.

```{r, eval = FALSE}
# Example: Analyze canopy height in Zurich using bounding box
res_bbox <- chm_analysis(
  bbox = c(8.51, 47.36, 8.56, 47.40),
  output_dir = tempdir(),
  max_tiles = 2,
  create_plots = TRUE
)
print(res_bbox$stats)

# Example: Use a location name (automatic geocoding)
res_loc <- chm_analysis(
  location = "Parc La Grange, Geneva, Switzerland",
  output_dir = tempdir(),
  max_tiles = 2
)
print(res_loc$stats)

# Example: Use your own canopy height raster
res_tif <- chm_analysis(
  chm_tif = "my_canopy.tif",
  output_dir = tempdir(),
  create_plots = TRUE,
  height_threshold = 3
)
print(res_tif$stats)
```

# 15. Urban Heat Island (UHI) Analysis

The `analyze_and_visualize_uhi()` function performs comprehensive Urban Heat Island analysis by integrating satellite-derived Land Surface Temperature (LST) data with environmental features from OpenStreetMap. This function automates the entire workflow from data acquisition to visualization.

**Key Features:**

- Automatic retrieval of thermal imagery (ECOSTRESS at 70m or Landsat 8/9 at 100m) via Microsoft Planetary Computer
- Integration with OSM data (green spaces, trees, buildings, water bodies)
- H3 hexagonal grid aggregation for consistent spatial analysis
- Getis-Ord Gi* hotspot analysis with statistical significance testing
- Moran's I spatial autocorrelation assessment
- Interactive Leaflet maps and publication-ready ggplot2 visualizations
- Export functions for GeoJSON, GeoPackage, CSV, and Shapefile formats

**Note:** This function requires internet access to download satellite imagery and OSM data.

```{r, eval = FALSE}
# Basic UHI analysis for a city
result <- analyze_and_visualize_uhi(
  location = "Basel, Switzerland",
  date_range = c("2023-06-01", "2023-08-31"),
  hex_resolution = 10,
  thermal_source = "auto"  # Tries ECOSTRESS first, then Landsat
)

# View the interactive map (opens in browser/viewer)
result$maps$interactive

# View static publication-ready maps
print(result$maps$combined)  # Combined 4-panel map
print(result$maps$scatter)   # Scatter plot dashboard

# Access statistics
print(result$stats$descriptive)
print(result$stats$correlations)

# Check metadata
print(result$meta)
```

The function returns a list containing:

- `results`: An sf object with hexagon-level metrics (LST, green coverage, built-up, hotspot categories)
- `maps`: Interactive Leaflet map and static ggplot2 maps
- `stats`: Descriptive statistics, correlations, regression results, and Moran's I
- `meta`: Processing metadata and data source information
- `export_geojson()` and `export_results()`: Functions for data export

```{r, eval = FALSE}
# Export results to multiple formats
result$export_results("basel_uhi", formats = c("geojson", "csv", "gpkg"))

# Using a bounding box instead of city name
result_bbox <- analyze_and_visualize_uhi(
  location = c(7.55, 47.53, 7.65, 47.58),  # xmin, ymin, xmax, ymax

  date_range = c("2023-07-01", "2023-07-31"),
  hex_resolution = 9  # Larger hexagons for faster processing
)

# With optional GHSL built-up data for more accurate built coverage
result_ghsl <- analyze_and_visualize_uhi(
  location = "Zurich, Switzerland",
  date_range = c("2023-06-01", "2023-08-31"),
  ghsl_path = "path/to/ghsl_built.tif",
  composite_scenes = TRUE,  # Composite multiple scenes to reduce cloud gaps
  max_scenes = 5
)
```

# 16. Shiny Application

You can make your own greenness analysis without having to code using an R Shiny implementation of the package. It is easily accessible from within R by calling the function `run_app()`.

# 17. Green View Index
This function allows the users to quantify urban greenness through image analysis. Utilizing the SuperpixelImageSegmentation library, it reads an image of an urban landscape and segments it into superpixels. The Green View Index (GVI) is then calculated by identifying green pixels within these segments. The GVI provides an objective measure of the proportion of visible vegetation in an image and is an important indicator for understanding urban greenness and its impact on ecological and human health.

```{r, eval = FALSE}
result <- calculate_and_visualize_GVI("image.png") #specify the image path
OpenImageR::imageShow(result$segmented_image) #To visualize the segmented image
green_pixels_raster <- as.raster(result$green_pixels_image) #To visualize green pixels
plot(green_pixels_raster)

# Save the segmented image
OpenImageR::writeImage(result$segmented_image, "segmented_image.png")

# Save the green pixels image
OpenImageR::writeImage(result$green_pixels_image, "green_pixels_image.png")
```

# 18. Green Space Similarity Index (GSSI)

The `gssi()` function calculates the Green Space Similarity Index (GSSI), a composite metric for evaluating and comparing urban green spaces across different regions by analyzing their size and spatial connectivity.
This function transforms spatial data into an equal-area projection for accurate area measurements, computes the total area of green spaces, and assesses their spatial connectivity using the Average Nearest Neighbor Distance (ANND). This dual approach provides a comprehensive view of the distribution and accessibility of green spaces.

The GSSI is calculated by inversely weighting the coefficient of variation in area sizes with the ANND, offering a score that reflects the abundance and accessibility of green spaces. Scores are normalized against the highest scoring city in the dataset for relative comparisons.


```{r, eval = FALSE}
d1 <- get_osm_data("New Delhi, India")
dsf <- d1$green_areas$osm_polygons
d2 <- get_osm_data("Basel, Switzerland")
bsf <- d2$green_areas$osm_polygons
d3 <- get_osm_data("Medellin, Colombia")
msf <- d3$green_areas$osm_polygons

cities_data <- list(dsf, bsf, msf)
gssi_values <- gssi(cities_data, "ESRI:54009")
```

# 19. Sky-View Factor (SVF) Microclimate Analysis

The `uh_svf()` function computes point-based and street-segment-level sky-view factors (SVF). SVF represents the horizontal projection of visible sky accounting for the Lambert cosine law of diffuse solar irradiance (Johnson & Watson 1984, Oke 1987). It calculates sky openness and obstruction from both terrain and buildings. 

It supports a hybrid workflow, automatically fetching building shapes from the Global Building Atlas S3 dataset, digital elevation models from AWS terrain tiles, and local canopy height models, or relying entirely on local datasets.

**Key Features:**
- Parallelized horizon-angle ray-casting across any number of azimuth directions.
- Generates interactive Leaflet maps with dynamic layers and linked legends.
- Generates fully interactive 3D WebGL explorers using Deck.gl extruded prisms.
- Generates publication-ready static cartographic maps and distribution curves.
- Computes radial horizon-line profiles to render Simulated Sky Obstruction Profiles for individual observer locations.

``` {r, eval = FALSE}
library(greenR)

# Basic SVF screening using AWS terrain and Global Building Atlas (GBA) building footprints
result_svf <- uh_svf(
  city_name = "Basel",
  bbox = c(7.58, 47.55, 7.60, 47.57),
  terrain_source = "elevatr",
  elevatr_z = 13,
  buildings_source = "gba",
  sample_mode = "street",
  spacing_street_m = 30,
  n_directions = 36,
  output_dir = tempdir(),
  output_prefix = "basel_svf"
)

# View the calculated microclimate statistics and summaries
print(result_svf$analytics)
```

# 20. Multi-Criteria Urban Heat and Canopy Decision Support Suite

The `uh_decision()` function is a high-level orchestration engine that automates the entire greenR Canyon Suite workflow. It compiles and merges multi-source spatial data to identify urban heat vulnerability and locate the highest-priority areas for greening interventions (such as tree-planting corridors).

It performs Multi-Criteria Decision Analysis (MCDA) based on Land Surface Temperature (LST), Population density, Sentinel-2 NDVI, Canopy Height Model (CHM) canopy cover, building footprints, and tree densities.

**Key Features:**
- Emulates street canyon orientations and solar bearings.
- Computes population-weighted spatial priority Gini inequality metrics with bootstrap confidence intervals.
- Segregates spatial priority zones into four actionable policy planning quadrants (e.g. "Hotter + Plantable").
- Generates 12 different publication-ready bivariate maps, 3D WebGL explorers, and multi-layer interactive Leaflet dashboards.

``` {r, eval = FALSE}
library(greenR)

# Run the complete decision support MCDA suite for a city
result_decision <- uh_decision(
  city_name = "Basel, Switzerland",
  hex_size_m = 100,
  ndvi_datetime = "2025-06-01/2025-08-31",
  lst_datetime = "2025-06-01/2025-08-31",
  output_dir = tempdir()
)

# Explore the resulting data layers
# Contains priority_grid (hex level), block_priority (superblock level), and canyon_priority (street corridor level)
print(names(result_decision))
```

