Skip to content

Commit 578cbf9

Browse files
committed
upmd
1 parent 27af520 commit 578cbf9

File tree

1 file changed

+142
-73
lines changed

1 file changed

+142
-73
lines changed

README.md

Lines changed: 142 additions & 73 deletions
Original file line numberDiff line numberDiff line change
@@ -15,18 +15,19 @@ This package is developed to facilitate robust and transparent analysis in green
1515
- [Installation](#installation)
1616
- [Functionalities](#functionalities)
1717
* [Availability](#availability)
18-
+ [Green Cover Streets](#green-cover-streets)
1918
+ [Calc NDVI](#calc-ndvi)
2019
+ [Land Cover](#land-cover)
2120
+ [Canopy coverage](#canopy-coverage)
2221
+ [Greenspace percentage](#greenspace-percentage)
22+
+ [Green Cover Streets](#green-cover-streets)
2323
* [Accessibility](#accessibility)
2424
+ [Greenspace access](#greenspace-access)
2525
* [Visibility](#visibility)
2626
+ [Viewshed](#viewshed)
2727
+ [vgvi from sf](#vgvi-from-sf)
2828
+ [vgvi from address](#vgvi-from-address)
29-
- [Extended Installation](#alternative-installation)
29+
- [Performance-errors-limitation-notes](#performance-errors-limitation-notes)
30+
- [Extended Installation](#extended-installation)
3031
* [GEE](#gee)
3132
* [Rcpp](#rcpp)
3233
- [Sources](#sources)
@@ -77,11 +78,11 @@ Please note that the examples based on this data serves as an illustration, and
7778

7879
Availability of greenness (reflecting the presence and amount of vegetation) will be assessed using four functions:
7980

80-
1. [Green Cover Streets](#green_cover_streets)
81-
2. [Calc NDVI](#calc-ndvi)
82-
3. [Land Cover](#land-cover)
83-
4. [Canopy coverage](#canopy-coverage)
84-
5. [Park percentage](#park-percentage)
81+
1. [Calc NDVI](#calc-ndvi)
82+
2. [Land Cover](#land-cover)
83+
3. [Canopy coverage](#canopy-coverage)
84+
4. [Park percentage](#park-percentage)
85+
5. [Green Cover Streets](#green_cover_streets)
8586

8687
Each function provide opportunity to measure greenness presence and amount in varying ways.
8788

@@ -90,37 +91,6 @@ The following subsections will briefly describe each availability function and e
9091

9192
---
9293

93-
### Green Cover Streets
94-
**[Global Coverage]** <br />
95-
The `green_cover_streets` function measure available greenery along streets in any city. Here we can use OSM street network for any city and connecting with European Space Agency's Worldcover map at 10m spatial resolution (https://esa-worldcover.org/en), we can measure how much trees, grass or shrubs are present around each street segment. We are going to use a buffer distance or a buffer zone of analysis around each street segment to estimate the percentage of tree, grass and shrubs and sum the total of these vegetation type to obtain total green coverage around each street segment. For that we are going to use the green_cover_streets() function.
96-
Here is an example for Amsterdam. To find the name of the place you are interested in try to use OSM address finder at exact place name using this link: https://nominatim.openstreetmap.org/ui/search.html
97-
98-
```r
99-
100-
#Measuring the green coverage around each street for the city, for a buffer distance of 30 m;
101-
# can try small area for fast implementation, try "Harmelen, Utrecht"
102-
103-
#for Cental Amsterdam
104-
GreenStcoverSt <- GreenExp::green_cover_streets ("Centrum, Amsterdam", buffer_distance = 30)
105-
106-
#map the result using the amazing mapview library for interactve mapping of the streets
107-
mapview::mapview(GreenStcoverSt, zcol = "greencover")
108-
109-
#Save the street file
110-
#To save the file we have to provide a path, for now, let us save it in current working diractory path
111-
path <- getwd() #change the path if you want to save it in specific folder
112-
113-
#save as in any GIS file format
114-
st_write(GreenStcoverSt, paste0(path,'/','GreenStcoverSt.gpkg'), delete_layer = TRUE)
115-
#can also be saved as shapefile, use GreenStcoverSt.shp
116-
117-
118-
```
119-
In the Figure below, we can find the Output after running the green_cover_streets function and a plot corresponding to the results.Here 0 values indicates no greenery within the buffer around the street, 100 indicates 100% presence of greenery within the given buffer zone around the street segment. We can change the background map to explore why some streets has high green coverage where others do not. Also on the interactive map we can click of each line to see what land cover type (i.e., tree, grass, shrub) contributes to overall green coverage. It also shows the presence of other land cover such as Cropland, which has not been considered for overall green coverage calculation.
120-
121-
<img src="man/figures/1_Green_streets.png" alt="Image" width="1000" />
122-
123-
---
12494

12595
### Calc NDVI
12696
**[Global Coverage]** <br />
@@ -313,6 +283,41 @@ The following figure shows the greenspace percentage value within 300m at each r
313283

314284
---
315285

286+
### Green Cover Streets
287+
**[Global Coverage]** <br />
288+
The `green_cover_streets` function measure available greenery along streets in any city. Here we can use OSM street network for any city and connecting with European Space Agency's Worldcover map at 10m spatial resolution (https://esa-worldcover.org/en), we can measure how much trees, grass or shrubs are present around each street segment. We are going to use a buffer distance or a buffer zone of analysis around each street segment to estimate the percentage of tree, grass and shrubs and sum the total of these vegetation type to obtain total green coverage around each street segment. For that we are going to use the green_cover_streets() function.
289+
Here is an example for Amsterdam. To find the name of the place you are interested in try to use OSM address finder at exact place name using this link: https://nominatim.openstreetmap.org/ui/search.html
290+
291+
```r
292+
293+
#Measuring the green coverage around each street for the city, for a buffer distance of 30 m;
294+
# can try small area for fast implementation, try "Harmelen, Utrecht"; "City of Melbourne", "Bogura, Bangladesh",
295+
# Large cities will take substential amount of time and computational power!
296+
297+
#for Cental Amsterdam
298+
GreenStcoverSt <- GreenExp::green_cover_streets ("Centrum, Amsterdam", buffer_distance = 30)
299+
300+
#map the result using the amazing mapview library for interactve mapping of the streets
301+
mapview::mapview(GreenStcoverSt, zcol = "greencover")
302+
303+
#Save the street file
304+
#To save the file we have to provide a path, for now, let us save it in current working diractory path
305+
path <- getwd() #change the path if you want to save it in specific folder
306+
307+
#save as in any GIS file format
308+
st_write(GreenStcoverSt, paste0(path,'/','GreenStcoverSt.gpkg'), delete_layer = TRUE)
309+
#can also be saved as shapefile, use GreenStcoverSt.shp
310+
311+
312+
```
313+
In the Figure below, we can find the Output after running the green_cover_streets function and a plot corresponding to the results.Here 0 values indicates no greenery within the buffer around the street, 100 indicates 100% presence of greenery within the given buffer zone around the street segment. We can change the background map to explore why some streets has high green coverage where others do not. Also on the interactive map we can click of each line to see what land cover type (i.e., tree, grass, shrub) contributes to overall green coverage. It also shows the presence of other land cover such as Cropland, which has not been considered for overall green coverage calculation.
314+
315+
<img src="man/figures/1_Green_streets.png" alt="Image" width="1000" />
316+
317+
---
318+
319+
---
320+
316321
## Accessibility
317322
The accessibility function of GreenExp is based on the concept of spatial proximity, meaning measuring distance to nearest greenspace.
318323
### Greenspace access
@@ -416,75 +421,114 @@ summary(building_access_Netj$greenspace_in_300m_buffer)
416421
#logical 2203 23437 251
417422

418423
```
419-
420-
424+
The results of these two accessibility approaches are considerably different. In particular, the network approach indicates longer median and mean distance compared to Euclidean distance. This is logical, considering the fact that street notwork distance are usually longer than crow-fly (Euclidean) distance. The functions provide different ways of measuring accessibility to ensure, access to greenspace has been estimated as with reasonable manner.
421425

422426

423427
## Visibility
424428
**[Local Coverage- Dependent on availability of input data]** <br />
429+
In GreenExp, we have functions to estimate visibility using Viewshed Greenness Visibility Index (VGVI) developed in [Labib et al., (2021)](https://doi.org/10.1016/j.scitotenv.2020.143050). The VGVI is a model that estimate potential greenness visibility using digital elevation data. Where presence of buildings or other obstacles can prevent the visible contact with greenery. The conceptual model of the VGVI index presented in figure below. The model details can be found at: [Labib et al., (2021)](https://doi.org/10.1016/j.scitotenv.2020.143050)
430+
<img src="man/figures/VGVI_model.jpg" alt="Image" width="1000" />
425431

426-
The visibility functions are made by the [GVI](https://github.com/STBrinkmann/GVI) package with some adaptations.
432+
Initially the VGVI was developed in python, later [Brinkmann & Labib, (2022)](https://github.com/STBrinkmann/GVI) created an R package [GVI](https://github.com/STBrinkmann/GVI) for fast calculation of the VGVI model. The visibility functions available in GreenExp package are from [GVI](https://github.com/STBrinkmann/GVI) package with some adaptations.
427433

428-
---
434+
**Note:** VGVI model needs LiDAR driven Digital Surface data, along with digital terrain data, which are usually unavailable for most places around the world. So the the model is only applicable in places where all three input data are available. The examples provided here are based on Dutch AHN dataset, which is freely available at national scale.
435+
436+
The fist function is simple viewshed analysis, which woks as a basis for VGVI modeling.
429437

430438
### Viewshed
431439

432-
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 usig 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
440+
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
433441
result of the `viewshed` function is a radial raster where 0 =
434442
no-visible and 1 = visible area.
443+
The code section below allow to use example data to run the viewshed and VGVI analysis
435444

436-
For a better explanation, go to the [GVI](https://github.com/STBrinkmann/GVI) package.
445+
``` r
437446

447+
#For this example, we first need to download the DSM, DTM and Green space layer from the following url (where we already uploaded the data)
438448

439-
```r
440-
# Read in the digital eleveation model
441-
GS <- terra::rast('data/GS_AMS.tif')
442-
# Read the digital surfaca model
443-
DSM <- terra::rast('data/DSM_AMS.tif')
444-
# Read the greenspace
445-
DTM<- terra::rast('data/DTM_AMS.tif')
449+
#create a temp diractory to get the zip file
450+
temp <- tempfile(fileext = "zip")
451+
download.file("https://a35b3ff7-2e06-4a4f-a669-8681322e59a7.usrfiles.com/archives/a35b3f_c9ec6cb349d043f4b3b055d6dc984ed3.zip", temp)
446452

447-
GreenExp::viewshed(observer = df_points[1,], dsm_rast = DSM, dtm_rast = DTM,
448-
max_distance = 200, observer_height = 1.7, plot = TRUE)
453+
#extract the zip file
454+
out <- unzip(temp, exdir = tempdir())
449455

450-
```
456+
#load the DSM, DTM, and GreenSpace (GS)
457+
DSM <- terra::rast(out[1])
458+
DTM <- terra::rast(out[2])
459+
GS <- terra::rast(out[3])
451460

452-
<img src="man/figures/viewshed.png" alt="Image" width="500" />
461+
#let us fix an obersver location
462+
obersver_mp <- sf::st_sf(sfheaders::sf_point(c(4.882752, 52.358029)), crs = st_crs(4326))
453463

454-
The left plot represents the Digital Elevation Model (DEM), whereas the right plot represents the viewshed, where green is the visible area and gray is not visible.
464+
#make sure the oberver location has same projection as the DSM/DTM/GS file. They all need to be in same projection
465+
obersver_mp <- sf::st_transform(obersver_mp, sf::st_crs(DSM))
455466

456-
### vgvi from sf
467+
singleview <- GreenExp::viewshed(observer = obersver_mp, dsm_rast = DSM, dtm_rast = DTM,
468+
max_distance = 200, observer_height = 1.7, plot = TRUE)
469+
470+
471+
```
472+
473+
474+
The result of viewshed model shows the visible and invisible cells from observer location.
457475

458-
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.
476+
<img src="man/figures/4_viewshed.png" alt="Image" width="500" />
459477

460-
Based on a viewshed and a binary greenspace raster, all visible points are classified as visible green and visible no-green. All values are summarised 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.
478+
The left plot represents the Digital Elevation Model (DEM), whereas the right plot represents the viewshed, where green is the visible area and gray is not visible.
461479

462-
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) and [Labib et al., 2021](https://doi.org/10.1016/j.scitotenv.2020.143050)
480+
### vgvi from sf
481+
Estimating visible greenery at point, line and polygon geometries using 'vgvi_from_sf' function. This function allows the user to provide multiple points, lines or polygons and then use those to create observer locations and estimate VGVI.
463482

464483
**Example:**
465484

466-
The figure below provides an overview of the data used for calculating the VGVI. In the figure, you will find various elements that contribute to the analysis.
485+
There are 3 main elements of VGVI modeling. These include:
467486

468487
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.
469488

470-
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.
489+
3. Green Space: The green here represent the a binary layer of greenery (in this case we are using a tree raster, 0= no tree, 1 = possibly a tree present in that raster cell), which represents the extent of green space. It indicates the areas covered by vegetation, such as trees or parks, and is a crucial factor in determining the VGVI.
471490

472-
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 parks, and is a crucial factor in determining the VGVI.
473-
474-
4. DEM and DSM: The figure also includes the Digital Elevation Model (DEM) and Digital Surface Model (DSM). These models provide information about the elevation of the terrain and structures present in the area. The combination of DEM and DSM helps in understanding the topography of the region.
491+
4. DEM and DSM: The Digital Elevation Model (DEM) and Digital Surface Model (DSM). These models provide information about the elevation of the terrain and structures present in the area. The combination of DEM and DSM helps in understanding the topography and viewing obstacles of the region.
475492

476493
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.
477494

478-
<img src="man/figures/vgvi_sf.png" alt="Image" width="500" />
495+
479496

480497
```r
481-
VGVI <- GreenExp::vgvi_from_sf(observer = df_points,
482-
dsm_rast = DSM, dtm_rast = DEM, greenspace_rast = GS,
498+
499+
# We need the same DSM, DTM and GS
500+
#If you have not downloaded them yet or they are not in R-environment, please load them first before running the code below
501+
502+
#for observer points, we can use street lines to generate sample points on streets. We can use OSM for that.
503+
# Here is an example for Central Amsterdam.
504+
505+
#get the bounding box
506+
boundbox <- tmaptools::geocode_OSM("Centrum, Amsterdam", as.sf = T, geometry = c("bbox"))
507+
508+
#get the lines
509+
lines <- osmextract::oe_get("Centrum, Amsterdam", stringsAsFactors=FALSE, boundary = boundbox,
510+
max_file_size = 5e+09, boundary_type = 'spat')
511+
512+
#Get the lines projected to input raster file, as they all have to be in same projection
513+
lines_prj <- sf::st_transform(lines, sf::st_crs(DSM))
514+
515+
516+
#Model VGVI on streets sample points, where spacing = 50 indicating the sample points will be at each 50 m
517+
VGVI_lines <- GreenExp::vgvi_from_sf(observer = lines_prj,
518+
dsm_rast = DSM, dtm_rast = DTM, greenspace_rast = GS, spacing = 50,
483519
max_distance = 200, observer_height = 1.7,
484-
m = 0.5, b = 8, mode = "logit")
520+
m = 0.5, b = 8, mode = "logit") #m and b are the distance decay parameter
521+
522+
#mapped results
523+
mapview::mapview(VGVI_lines, zcol = "VGVI")
524+
525+
#VGVI points over lines
526+
mapview::mapview(VGVI_lines, zcol = "VGVI") + mapview::mapview(lines_prj)
527+
485528
```
529+
The result of the above function shows below. Here 0 indicates no tree visibility, 1 indicates 100% tree visibility within 200m viewing distance from each over location. Note these these points can be also joined back to street again for obtaining aggregated tree visibility on each road segment.
486530

487-
531+
<img src="man/figures/4_VGVI_lines.png" alt="Image" />
488532

489533
---
490534

@@ -495,19 +539,44 @@ GVI from address function: In contrast, the VGVI from address function employs a
495539

496540
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 GVI within a defined buffer and provides a more averaged assessment of the VGVI.
497541

542+
For an example for some random locations in Amsterdam we can estimate the VGVI based on samples around them
543+
498544

499545
```r
500-
VGVI <- GreenExp::vgvi_from_address(address = df_points,
501-
dsm_rast = DSM, dtm_rast = DEM, greenspace_rast = GS,
546+
547+
#get the OSM city geocoded bounding box
548+
getcityboudingbox <- tmaptools::geocode_OSM("Centrum, Amsterdam", as.sf = T, geometry = c("bbox"))
549+
550+
#generate random points within the bounding box
551+
RandomObservePoints <- sf::st_sample(getcityboudingbox, size = 1000) %>% st_as_sf()
552+
553+
#project the observer points to match with the DSM/DTM/GS
554+
RandomObservePoints <- sf::st_transform(RandomObservePoints, sf::st_crs(DSM))
555+
556+
557+
VGVI_random_points <- GreenExp::vgvi_from_address(address = RandomObservePoints,
558+
dsm_rast = DSM, dtm_rast = DTM, greenspace_rast = GS, buffer_distance = 50,
502559
max_distance = 200, observer_height = 1.7,
503560
m = 0.5, b = 8, mode = "logit")
504561

562+
#map the results
563+
mapview::mapview(VGVI_random_points, zcol = "mean_VGVI")
564+
565+
505566
```
567+
The results showing mean VGVI value for each points based on sample points generated around them.
568+
<img src="man/figures/4_VGVI_randompoints.png" alt="Image" />
506569

507-
<img src="man/figures/vgvi_address.png" alt="Image" width="500" />
508570

509571
---
510572

573+
# Performance-errors-limitation-notes
574+
GreenExp is in still development phase, and have ongoing issues that we are progressing. So the package has some performance issues and limitations. Here are some common limitations and errors might happen when you are testing the package. You can help improving it in future!
575+
- Running Speed: Some large areas extracted from OSM, might take a long time to run different metrics. Such as the street green cover function takes a long time for very large area. We suggest you can test it in small areas.
576+
- Missing
577+
578+
579+
511580
# Extended installation
512581

513582
To make optimal use of the package

0 commit comments

Comments
 (0)