| Type: | Package | 
| Title: | Access and Analyze Global GreenSpace Spatial Data | 
| Version: | 0.1.1 | 
| Description: | Access and analyze multi-band greenspace seasonality data cubes (available for 1,028 major global cities), global Normalized Difference Vegetation Index / land cover data from the European Space Agency WorldCover 10m Dataset, and Sentinel-2-l2a images. Users can download data using bounding boxes, city names, and filter by year or seasonal time window. The package also supports calculating human exposure to greenspace using a population-weighted greenspace exposure model introduced by Chen et al. (2022) <doi:10.1038/s41467-022-32258-4> based on Global Human Settlement Layer population data, and calculating a set of greenspace morphology metrics at patch and landscape levels. | 
| URL: | https://github.com/billbillbilly/greenSD, https://billbillbilly.github.io/greenSD/ | 
| BugReports: | https://github.com/billbillbilly/greenSD/issues | 
| License: | GPL-3 | 
| Encoding: | UTF-8 | 
| Depends: | R (≥ 4.1) | 
| Suggests: | testthat (≥ 3.0.0), knitr, rmarkdown | 
| VignetteBuilder: | knitr | 
| Imports: | terra, maptiles, sf, future, furrr, purrr, nominatimlite, utils, cli, dsmSearch, rstac, landscapemetrics, magick, aws.s3, dplyr, tidyr, stringr, tibble, grDevices, rlang | 
| RoxygenNote: | 7.3.2 | 
| NeedsCompilation: | no | 
| Packaged: | 2025-10-26 23:24:16 UTC; xiaohaoyang | 
| Author: | Xiaohao Yang [aut, cre, cph] | 
| Maintainer: | Xiaohao Yang <xiaohaoy111@gmail.com> | 
| Repository: | CRAN | 
| Date/Publication: | 2025-10-30 19:50:02 UTC | 
Get all of the urban areas in the Greenspace Seasonality Data Cube
Description
This function returns all of the urban areas in the Greenspace Seasonality Data Cube dataset.
Usage
check_available_urban(test = FALSE)
Arguments
| test | logical. (ignored) Only for testing. | 
Value
dataframe
Note
You can explore all available urban areas in an interacive map at: https://github.com/billbillbilly/greenSD/blob/main/scripts/city_urban_boundaries.geojson
References
Wu, S., Song, Y., An, J. et al. High-resolution greenspace dynamic data cube from Sentinel-2 satellites over 1028 global major cities. Sci Data 11, 909 (2024). https://doi.org/10.1038/s41597-024-03746-7
Examples
check_available_urban(test = TRUE)
Get an urban area boundary based on the UID
Description
This function returns a polygon of a city boundary based on the UID
Usage
check_urban_boundary(uid = NULL, plot = TRUE, test = FALSE)
Arguments
| uid | numeric. Urban area ID. To check the ID of an available urban area,
use  | 
| plot | logical. Whether to plot city boundary | 
| test | logical. (ignored) Only for testing. | 
Value
sf
References
Wu, S., Song, Y., An, J. et al. High-resolution greenspace dynamic data cube from Sentinel-2 satellites over 1028 global major cities. Sci Data 11, 909 (2024). https://doi.org/10.1038/s41597-024-03746-7
Examples
check_urban_boundary(test = TRUE)
exposure
Description
Computes population-weighted greenspace fraction or human exposure to greenspace based on a population-weighted exposure model (Chen et al., 2022), using population data from the Global Human Settlement Layer (GHSL; Pesaresi et al., 2024). See Details for the underlying method and assumptions.
Usage
compute_exposure(
  r = NULL,
  res = c(10, 10),
  pop_year = 2020,
  radius = 500,
  grid_size = NULL,
  height = FALSE,
  pop_out = FALSE,
  quiet = TRUE
)
Arguments
| r | A SpatRaster with single/multiple greenspace layer(s), either
fractional or binary (where non-green = 0 and green = 1), typically
the output from  | 
| res | numeric vector of length 2. The actual spatial resolution (in meters).
Default is  | 
| pop_year | numeric. Year of the GHSL dataset to use. Must be one of: 2015, 2020, 2025, or 2030. Default is 2020. | 
| radius | numeric. Buffer radius (in meters) used for local averaging.
Default is  | 
| grid_size | numeric. Optional. If provided, output is aggregated to grid cells
of this size (in meters) and returned as an  | 
| height | logical. Whether to compute greenspace volume for population-weighted greenspace fraction or human exposure to greenspace using Meta's global canopy height map (Tolan et al., 2024). (The default is FALSE) | 
| pop_out | logical. Whether return population layer. | 
| quiet | logical. Whether show progress bars for some process. | 
Details
This function implements the population-weighted greenspace exposure (PWGE) model:
- Start with a population raster. Each pixel \( i \) has a population value \( P_i \). 
- Create a circular buffer of radius \( d \) around each pixel center. 
- For each buffer, calculate greenspace fraction: - G_i^d = \frac{\text{Area of greenspace within buffer}}{\text{Total buffer area}}
- Repeat for all \( i = 1, 2, ..., N \) grid cells. 
- Compute overall exposure: - GE^d = \frac{\sum_i P_i \cdot G_i^d}{\sum_i P_i}
Value
SpatRaster or sf. A SpatRaster (if grid_size is NULL) with
layers pwgf_*, or an sf object with columns pwge_* representing
population-weighted greenspace exposure values aggregated to each grid polygon.
References
Chen, B., Wu, S., Song, Y. et al. Contrasting inequality in human exposure to greenspace between cities of Global North and Global South. Nat Commun 13, 4636 (2022). https://doi.org/10.1038/s41467-022-32258-4
Pesaresi, M., Schiavina, M., Politis, P., Freire, S., Krasnodębska, K., Uhl, J. H., … Kemper, T. (2024). Advances on the Global Human Settlement Layer by joint assessment of Earth Observation and population survey data. International Journal of Digital Earth, 17(1). https://doi.org/10.1080/17538947.2024.2390454
Tolan, J., Yang, H. I., Nosarzewski, B., Couairon, G., Vo, H. V., Brandt, J., ... & Couprie, C. (2024). Very high resolution canopy height maps from RGB imagery using self-supervised vision transformer and convolutional decoder trained on aerial lidar. Remote Sensing of Environment, 300, 113888.
Examples
sample_data <- terra::rast(system.file("extdata", "detroit_gs.tif", package = "greenSD"))
pwgf <- compute_exposure(
  # r = sample_data,
  pop_year = 2020,
  radius = 1500
)
compute_morphology
Description
Compute greenspace morphology metrics at patch (Nowosad & Stepinski, 2019) or landscape level (see details), including average size (AREA_MN), fragmentation (PD), connectedness (COHESION), aggregation (AI), and complexity of the shape (SHAPE_AM), related to public health (Wang et al., 2024)
Usage
compute_morphology(r = NULL, directions = 4, grid_size = NULL, quiet = TRUE)
Arguments
| r | SpatRaster. A single-band binary greenspace raster, where 0 or NA represents non-green areas and 1 represents green areas. | 
| directions | numeric. The number of directions in which patches should be connected: 4 (default) or 8. | 
| grid_size | numeric or sf polygons. (Optional) If specified, morphology metrics at grid level will be computed based on the size (in meters) of given grid cells or input (sf) polygons. | 
| quiet | logical. Whether show progress bars for some process. | 
Details
To get information of metrics, please use landscapemetrics::list_lsm().
Value
A SpatVector object contains indivisual patches with metrics at patch level,
when grid_size = NULL.
A SpatVector object contains landscape-level value of metrics,
when grid_size is not NULL.
References
Nowosad J., TF Stepinski. 2019. Information theory as a consistent framework for quantification and classification of landscape patterns. https://doi.org/10.1007/s10980-019-00830-x
Wang, H., & Tassinary, L. G. (2024). Association between greenspace morphology and prevalence of non-communicable diseases mediated by air pollution and physical activity. Landscape and Urban Planning, 242, 104934.
Examples
green <- get_tile_green(
                        # bbox = c(-83.087174,42.333373,-83.042542,42.358748),
                        provider = "esri",
                        zoom = 16)
# p <- terra::ifel(green$green == 0, NA, 1)
m <- compute_morphology(
                       #r = p
                       directions = 8)
Get band index based on time period
Description
Converts a date string in "MM-DD" format to the corresponding band index
for the Greenspace Seasonality Data Cube, which contains 36 bands representing
10-day intervals over a year.
Usage
get_band_index_by_time(time, year)
Arguments
| time | Character vector of length 2. (optional) Start and end dates in  | 
| year | numeric. (required) The year of interest. | 
Details
The Greenspace Data Cube is organized into 36 bands per year, each representing a 10-day interval. This function calculates which of those bands a given date falls into by converting the MM-DD string into the day-of-year (DOY) and dividing by 10 (rounded up).
Value
Interger. a band index.
Examples
get_band_index_by_time(c("03-20", "10-15"), year = 2020)
Download landcover or NDVI Data from ESA WorldCover 10m Annual Dataset
Description
download 11-class landcover or 3-band NDVI Data (NDVI p90, NDVI p50, NDVI p10). Users can define an area of interest using a bounding box or place name.
Usage
get_esa_wc(
  bbox = NULL,
  place = NULL,
  datatype = "landcover",
  year = 2021,
  mask = TRUE,
  quiet = TRUE
)
Arguments
| bbox | 
 | 
| place | character or vector. (optional) A single line address, e.g. ("1600 Pennsylvania Ave NW, Washington") or a vector of addresses (c("Madrid", "Barcelona")). | 
| datatype | character. One of "landcover" and "ndvi". | 
| year | numeric. The year of interest:  | 
| mask | logical (optional). Default is  | 
| quiet | logical. Whether show progress bars for some process. | 
Value
A SpatRaster object containing 11-class land cover or NDVI
yearly percentiles composite (NDVI p90, NDVI p50, NDVI p10)
References
Zanaga, D., Van De Kerchove, R., De Keersmaecker, W., Souverijns, N., Brockmann, C., Quast, R., Wevers, J., Grosu, A., Paccini, A., Vergnaud, S., Cartus, O., Santoro, M., Fritz, S., Georgieva, I., Lesiv, M., Carter, S., Herold, M., Li, L., Tsendbazar, N.-E., … Arino, O. (2021). ESA WorldCover 10 m 2020 v100 (Version v100). Zenodo. https://doi.org/10.5281/zenodo.5571936
Zanaga, D., Van De Kerchove, R., Daems, D., De Keersmaecker, W., Brockmann, C., Kirches, G., Wevers, J., Cartus, O., Santoro, M., Fritz, S., Lesiv, M., Herold, M., Tsendbazar, N.-E., Xu, P., Ramoino, F., & Arino, O. (2022). ESA WorldCover 10 m 2021 v200 (Version v200). Zenodo. https://doi.org/10.5281/zenodo.7254221
Examples
result <- get_esa_wc(
  # place = 'New York'
  year = 2021
)
Download Greenspace Seasonality Data Cube
Description
download Greenspace Seasonality Data Cube for an urban area. Retrieves high-resolution greenspace seasonality data from the Sentinel-2-based global dataset developed by Wu et al. (2024). Users can define a city of interest using a bounding box, place name, coordinates, or unique city ID (UID).
Usage
get_gsdc(
  bbox = NULL,
  place = NULL,
  location = NULL,
  UID = NULL,
  year = NULL,
  time = NULL,
  mask = TRUE,
  quiet = TRUE
)
Arguments
| bbox | 
 | 
| place | character or vector. (optional) A single line address,
e.g. ("1600 Pennsylvania Ave NW, Washington") or a vector of addresses
(c("Madrid", "Barcelona")). This can be ignored if  | 
| location | vector or sf point. A point of interest.
Ignored if  | 
| UID | numeric. Urban area ID. To check the ID of an available urban area,
use  | 
| year | numeric. (required) The year of interest. | 
| time | Character vector of length 2 or character. (optional) Start and end dates in
 | 
| mask | logical (optional). Default is  | 
| quiet | logical. Whether show progress bars for some process. | 
Details
The Greenspace Data Cube is organized into 36 bands per year, each representing a 10-day interval.
Value
A SpatRaster object containing the greenspace seasonality data.
Note
Use check_available_urban() and check_urban_boundary() to see supported
cities and their boundaries.
References
Wu, S., Song, Y., An, J. et al. High-resolution greenspace dynamic data cube from Sentinel-2 satellites over 1028 global major cities. Sci Data 11, 909 (2024). https://doi.org/10.1038/s41597-024-03746-7
Examples
result <- get_gsdc(UID = 0,
                   # year = 2022
                  )
Retrieve Sentinel-2-l2a images to compute NDVI
Description
download Sentinel-2-l2a imagery data and compute NDVI. Users can define an area of interest using a bounding box or place name.
Usage
get_s2a_ndvi(
  bbox = NULL,
  place = NULL,
  datetime = c(),
  cloud_cover = 10,
  vege_perc = 0,
  select = "latest",
  method = "first",
  mask = TRUE,
  output_bands = NULL,
  quiet = TRUE
)
Arguments
| bbox | 
 | 
| place | character or vector. (optional) A single line address, e.g. ("1600 Pennsylvania Ave NW, Washington") or a vector of addresses (c("Madrid", "Barcelona")). | 
| datetime | numeric vector of 2. The time of interest such as
 | 
| cloud_cover | numeric. Threshold for the percentage of cloud coverage. Desfault is 10. | 
| vege_perc | numeric. Threshold for the percentage of vegetation coverage. Desfault is 0. | 
| select | character. one of "latest", "earliest", "all". The default is "latest". | 
| method | character. A method for mosaicing layers: one of "mean", "median", "min", "max", "modal", "sum", "first", "last". The default is "first". | 
| mask | logical (optional). Default is  | 
| output_bands | vector. A list of band names ( | 
| quiet | logical. Whether show progress bars for some process. | 
Value
A SpatRaster object containing (multiple) NDVI layer(s) (for different
period of time) select = "latest" or select = "first"
(or if mask = TRUE and select = "all")
A List of NDVI rasters if mask = FALSE and select = "all".
Examples
result <- get_s2a_ndvi(
  # place = 'New York',
  datetime = c("2020-08-01", "2020-09-01")
)
Classify greenspace based on map tile images
Description
Generate high-resolution greenspace segmentation using WorldImagery map tiles provided by esri and Sentinel-2 cloudless mosaic tiles provided by EOX.
Usage
get_tile_green(
  bbox = NULL,
  place = NULL,
  zoom = 17,
  provider = "esri",
  year = NULL,
  quiet = TRUE
)
Arguments
| bbox | 
 | 
| place | character or vector. (optional) A single line address, e.g. ("1600 Pennsylvania Ave NW, Washington") or a vector of addresses (c("Madrid", "Barcelona")). | 
| zoom | numeric. Zoom level of map tile. The default is  | 
| provider | character. One of "esri" and "eox". | 
| year | integer. The desired year for Sentinel-2 cloudless mosaic
tiles. (This is required when  | 
| quiet | logical. Whether show progress bars for some process. | 
Value
A list of two rasters including: greenspace segmentation (where 1 is green and 0 is non-green) and original map tiles
Note
The data derived from Esri WorldImagery may need to include appropriate Esri copyright notice.
Examples
g <- get_tile_green(
 # bbox = c(-83.087174,42.333373,-83.042542,42.358748),
 zoom = 15
)
ndvi_to_sem
Description
Convert ndvi raster data into semantic vegetation areas
Usage
ndvi_to_sem(r = NULL, threshold = c(0.2, 0.5), quiet = FALSE)
Arguments
| r | A SpatRaster with single greenspace layer, typically
the output from  | 
| threshold | numeric vector of two. Thresholds, defaulting to  | 
| quiet | logical. Whether show progress bars for some process. | 
Value
SpatRaster. A raster, where 0 represents non-green area, 1 represents
shrub and grassland, and 2 represents trees.
References
Hashim, H., Abd Latif, Z., & Adnan, N. A. (2019). Urban vegetation classification with NDVI threshold value method with very high resolution (VHR) Pleiades imagery. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 42, 237-240.
Examples
sample_data <- terra::rast(system.file("extdata", "detroit_gs.tif", package = "greenSD"))
seg <- ndvi_to_sem(sample_data$`25_NDVI`, threshold = c(0.2, 0.6))
Sample greenspace-realted data from Greenspace Seasonality Data Cube, ESA WorldCover 10m Annual Composites Dataset, or Sentinel-2-l2a images.
Description
Samples values by locatoins from the Greenspace Seasonality Data Cube developed by Wu et al. (2024), ESA WorldCover 10m Annual Composites Dataset by Zanaga et al. (2021), or Sentinel-2-l2a images.
Usage
sample_values(
  samples = NULL,
  time = NULL,
  source = "gsdc",
  output_bands = NULL,
  cloud_cover = 10,
  vege_perc = 0,
  select = "latest",
  method = "first",
  quiet = TRUE
)
Arguments
| samples | A list, matrix,  | 
| time | numeric or vector. The time of interest. See Detail. | 
| source | character. The data source for extracting greenspace values:
 | 
| output_bands | vector. A list of band names ( | 
| cloud_cover | numeric. The percentage of cloud coverage for retrieving
Sentinel-2-l2a images. (Only required, when  | 
| vege_perc | numeric. The percentage of cloud coverage for retrieving
Sentinel-2-l2a images. (Only required, when  | 
| select | character. one of "latest", "earliest", "all". The default is "latest". | 
| method | character. A method for mosaicing layers: one of "mean", "median", "min", "max", "modal", "sum", "first", "last". The default is "first". | 
| quiet | logical. Whether show progress bars for some process. | 
Details
time: For the greenspace seasonality data cube, only years from 2019 to 2022
are availabe. For ESA WorldCover 10m Annual Composites Dataset, only 2020
and 2021 are available.
Value
A data.frame containing greenspace values extracted at each point
across all bands. Each row corresponds to a sample location;
columns represent band values.
Note
For sampling data from Greenspace Seasonality Data Cube samples must be
located within the same boundary of an available city in the data cube.
Use check_available_urban() and check_urban_boundary() to see supported
cities and their boundaries.
References
Wu, S., Song, Y., An, J. et al. High-resolution greenspace dynamic data cube from Sentinel-2 satellites over 1028 global major cities. Sci Data 11, 909 (2024). https://doi.org/10.1038/s41597-024-03746-7
Zanaga, D., Van De Kerchove, R., De Keersmaecker, W., Souverijns, N., Brockmann, C., Quast, R., Wevers, J., Grosu, A., Paccini, A., Vergnaud, S., Cartus, O., Santoro, M., Fritz, S., Georgieva, I., Lesiv, M., Carter, S., Herold, M., Li, L., Tsendbazar, N.-E., … Arino, O. (2021). ESA WorldCover 10 m 2020 v100 (Version v100). Zenodo. https://doi.org/10.5281/zenodo.5571936
Zanaga, D., Van De Kerchove, R., Daems, D., De Keersmaecker, W., Brockmann, C., Kirches, G., Wevers, J., Cartus, O., Santoro, M., Fritz, S., Lesiv, M., Herold, M., Tsendbazar, N.-E., Xu, P., Ramoino, F., & Arino, O. (2022). ESA WorldCover 10 m 2021 v200 (Version v200). Zenodo. https://doi.org/10.5281/zenodo.7254221
Examples
# see supported urban areas and their boundaries
check_available_urban()
boundary <- check_urban_boundary(uid = 11)
# sample locations with in the boundary
samples <- sf::st_sample(boundary, size = 20)
# extract values
gs_samples <- sample_values(samples,
                            # time = 2022
                           )
Convert A Multi-layer Raster to GIF
Description
Export a multi-layer raster (SpatRaster) or vector layer (sf)
with multiple numeric value columns to an animated GIF.
Usage
to_gif(
  r,
  fps = 5,
  width = 600,
  height = 600,
  axes = TRUE,
  title_prefix = NULL,
  border = FALSE
)
Arguments
| r | SpatRaster or sf. A SpatRaster with multiple layers or an sf object with multiple numeric value columns. | 
| fps | numeric. Frames per second (default 5). | 
| width | numeric. Width of output GIF in pixels. | 
| height | numeric. Height of output GIF in pixels. | 
| axes | logical. Draw axes? | 
| title_prefix | character or character vector. | 
| border | character. Color of polygon border(s); using NA hides them.
Only optional when  | 
Value
An animated magick image object (GIF).
Examples
sample_data <- terra::rast(system.file("extdata", "detroit_gs.tif", package = "greenSD"))
gif <- to_gif(sample_data)