Tree Data International Comparison Documentation

International Comparison Documentation

This page attempts to show how the international comparison was performed, using Vienna as an example.

The first step is to get a 20km by 20km radius around the city of choice, for logistical reasons this is a square.

Code

#build dataframes

library(stars)
library(raster)
library(sf)
library(dplyr)
library(leaflet)

coordinates <- list(c(16.3713,48.2081 )) 

crss <- c(31287)

cities <- c( 'vienna')
       

source('../r/functions.r')
source('../r/mapping_functions.r')
source('../themes/theme.r')

#settings
radius <- 10000
union <- TRUE
cell_size <- 1000
res <- 5
threshold <- 0

#start actual code
city_df = tibble(coords = coordinates, coordinateref = crss, city = cities) %>%
  as.data.frame() %>%
  rowwise() %>%
  mutate(lng = coords[1], lat = coords[2])

# city_df = city_df[-c(1:(nrow(city_df) -3 ) ),]
# 
# city_df = city_df %>% filter(city == 'osaka')

# city_df = city_df %>%
#   rowwise() %>%
#   mutate(res = runCity(lng, lat, coordinateref, city) )

#get water layer ready

water_bodies_global <- read_sf('../data/water_valid.shp')

  longit = city_df$lng[1]
  latit = city_df$lat[1]
  local_crs = city_df$coordinateref[1]
  cur_city = city_df$city[1]

  local_coords = c(as.numeric(longit),as.numeric(latit))
  
  centre = st_point(local_coords) %>%
    st_sfc(crs = 4326) 
  
  sf_use_s2(F)

  bounding = centre %>%
    st_transform(crs = local_crs) %>%
    st_buffer(dist = radius) %>%
    st_bbox() %>%
    st_as_sfc(crs = 4326)
  
  bdmap <- bounding %>%
    st_transform('wgs84') %>%
    leaflet() %>%
    addProviderTiles('CartoDB.Positron') %>%
    addMeasure(primaryLengthUnit = 'meters') %>%
    addPolygons()
  
  
  bdmap

We then need to find where the water is as to not disrupt the estimates with zeros.

Code

  water_bodies = st_transform(water_bodies_global, local_crs)

  #water_bodies = water_bodies %>% filter(TYPE == 'Ocean or Sea')

  wbm = st_is_valid(water_bodies)
  wb = water_bodies %>% filter( wbm )

  sf_use_s2(F)

  ww <- st_intersection(wb, bounding)
  
  ww <- st_as_sf(ww) %>%
    dplyr::select(SHAPE_Area)
  
  ww = distinct(ww) %>%
    rowwise() %>%
    mutate(area = st_area(geometry)) %>%
    ungroup() %>%
    mutate(bounding_area = st_area(bounding)) %>%
    mutate(SHAPE_Area = round(SHAPE_Area, 5)) %>%
    rowwise() %>%
    mutate(isNotEqlToBounding = !(bounding_area == area )) %>%
    filter(isNotEqlToBounding) %>%
    ungroup() %>%
    filter(SHAPE_Area != 0.06180, SHAPE_Area != 0.05493)
  
  wwu = st_union(ww)

  water_map <- wwu %>%
    st_transform('wgs84') %>%
    leaflet() %>%
    addProviderTiles('CartoDB.Positron') %>%
    addPolygons()
  
  water_map

After finding where the water is, we need to then cut up the city into a grid, each 1km by 1km. This gives us 400 ‘chunks’ to perform the analysis. Now we take a look at an example grid square.

Code

  grid <- st_make_grid(bounding, cellsize = c(cell_size,cell_size), square = TRUE)
  n_cells <- grid %>% length()
  
  
  grid %>%
    st_transform(crs = "wgs84") %>%
    leaflet() %>%
    addMeasure(primaryLengthUnit = 'meters') %>%
    addProviderTiles('CartoDB.Positron') %>%
    addPolygons()

First, find the bonding box for each grid square.

Code

  index = 272
  bbox = grid[index,]
  city = 'vienna'
  water = wwu
   
  if( !dir.exists( paste0('../rasters/international/', city)) ) {
    dir.create( paste0('../rasters/international/', city) )
  }
  
  
  distance_to_centre <- st_distance( centre, st_centroid( st_transform( bbox , 4326) ) )
  
  item_bbox <- bbox %>%
    st_transform(4326)
  
  mp(item_bbox)

Then we subtract the water layer.

Code

  if( !dir.exists( paste0('../rasters/international/',city, '/',  'rad', radius, 'cellsize',cell_size, res,'res/')) ) {
    dir.create( paste0('../rasters/international/',city, '/',  'rad', radius, 'cellsize',cell_size, res,'res/') )
  }
  
  file_name = paste0('../rasters/international/',city, '/',  'rad', radius, 'cellsize',cell_size, res,'res/',index,'.tif')
  
  #print(file_name)
  
  if(file.exists(file_name)) {
    #print('running from file')
    
    chm <- raster(file_name)
    var = paste0('X',index)
    
  } else {
    print('downloading from api')
    chm <- chmloader::download_chm(item_bbox, filename = file_name, res = res)
    var = paste0(index,'.tif')
  }
  
  
  chm_stars = st_as_stars(chm, proxy = T )
  chm_simplified = st_as_sf(chm_stars, merge = TRUE) 
  
  chm_simplified = chm_simplified %>% filter( !!as.name(var) > threshold) 
  chm_simplified <- st_transform(st_as_sf(chm_simplified), 7844)
  
  chm_is_empty = ifelse(nrow(chm_simplified) > 0, F, T) 
  
  if(!chm_is_empty) {
    if(union && nrow(chm_simplified) > 0){
      chm_simplified = st_union(chm_simplified)
    }
  }
  
  bbox_no_water = st_difference(wgs(bbox), wgs(wwu))
  
  mp(wgs(bbox_no_water))

Then, we find the tree coverage for that square, excluding water. We can then calculate a simple percentage covering of the grid square.

Code

  chm_simplified = st_intersection(wgs(chm_simplified), wgs(bbox_no_water))

although coordinates are longitude/latitude, st_intersection assumes that they
are planar

Code

  mp(wgs(chm_simplified))

Code

  no_water_area = st_area(bbox_no_water) %>%
    as.vector()
  
  pc <- st_area(chm_simplified) / no_water_area
  
  pc = round(pc, 4)

Then we get our final answer of 0.3154. We then repeat for all 400 grid items, and then save the results to a file for later analysis.