Spatial Data Analysis (450:320)

Class 13

Today

  • Attribute operations
  • Geometric measurements
  • Geometric operations
  • Geometry operations

Use osmdata to query some vectors

library(dplyr)
library(osmdata) # has to install first

q <- opq(bbox = st_bbox(
    c(xmin = -74.448598, xmax = -74.429830, 
      ymax = 40.532277, ymin = 40.518258), crs = st_crs(4326)))

# Building
buildings <- q %>% 
  add_osm_feature(key = 'building') %>%
  osmdata_sf() %>% `[[`('osm_polygons') %>% 
  filter(!is.na(name))

# Streets
streets <- q %>% 
  add_osm_feature(key = 'highway') %>%
  osmdata_sf() %>% `[[`('osm_lines') %>% 
  filter(!is.na(name))

# Trees
trees <- q %>% 
  add_osm_feature(key = c('natural')) %>%
  osmdata_sf() %>% `[[`('osm_points')

The data

plot(buildings$geometry, col = "black")
plot(streets %>% st_geometry(), add = T, col = "brown")
plot(st_geometry(trees), pch = 20, add = TRUE, col = "darkgreen")

Attribute operations

  • sf is a data.frame, so all base R and dplyr syntax works on it. E.g.

    • subsetting
    • aggregation (group)
    • joining

  • st_drop_geometry() to drop geometry column

Practice 1

Tasks:

  • Select out streets that are residential or path.
  • Count the number of streets by types.
# A tibble: 6 × 2
  highway      count
  <chr>        <int>
1 path             8
2 residential     50
3 secondary       11
4 service          2
5 tertiary        19
6 unclassified     1

Geometric measurements

  • st_area
  • st_length
  • st_distance

st_area

buildings <- buildings %>% 
    mutate(area = st_area(.)) %>% 
    select(name, area)

head(buildings, 3) # Observe column area
Simple feature collection with 3 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -74.43991 ymin: 40.5225 xmax: -74.4353 ymax: 40.52598
Geodetic CRS:  WGS 84
                             name           area                       geometry
221488774 Livingston Apartments C 4513.754 [m^2] POLYGON ((-74.4367 40.52485...
221488775 Livingston Apartments B 3942.836 [m^2] POLYGON ((-74.43778 40.5256...
251381247               Beck Hall 1841.321 [m^2] POLYGON ((-74.43915 40.5225...
ggplot(buildings) +
    geom_sf(color = "black", aes(alpha = as.numeric(area))) +
    labs(alpha = "Area (m^2)") +
    theme_void()

Units

units::set_units()

# Examples of usage
x <- 10 %>% units::set_units('km2')
x
10 [km^2]
# Remove units
x <- units::set_units(x, NULL)
x
[1] 10
buildings$area <- units::set_units(buildings$area, 'ha')
head(buildings, 3)
Simple feature collection with 3 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -74.43991 ymin: 40.5225 xmax: -74.4353 ymax: 40.52598
Geodetic CRS:  WGS 84
                             name           area                       geometry
221488774 Livingston Apartments C 0.4513754 [ha] POLYGON ((-74.4367 40.52485...
221488775 Livingston Apartments B 0.3942836 [ha] POLYGON ((-74.43778 40.5256...
251381247               Beck Hall 0.1841321 [ha] POLYGON ((-74.43915 40.5225...

Practice 2

streets:

  • Calculate the length of each street in km and save to a new column “length” to streets.
  • Calculate the average length by street types (in column highway).

trees: Calculate the distance between the first two tree in miles.

# A tibble: 6 × 2
  highway      avg_length
  <chr>              [km]
1 path              0.649
2 residential       0.285
3 secondary         0.122
4 service           0.165
5 tertiary          0.157
6 unclassified      0.486

Geometric operations

  • st_centroid (get the centroid of objects)
  • st_buffer (add buffer to any objects)
  • st_coordinates (get coordinates of objects)
  • st_simplify (reduce the complexity of objects)
  • st_bbox (get the bounding box of objects)

Practice 3

  • Simplify the buildings and get the centroid of the them.
  • Reproject streets to Albers Equal Area (EPSG:5070) and add a 10 meter buffer to them.
  • Get the coordinates of trees.

Geometry operations

  • Clipping/overlapping

    • st_intersection
    • st_difference & st_sym_difference
    • st_union: return sfc

Practice 4

# Your building from osmdata
building_osm <- buildings %>% filter(name == '...')

# building polygon you created
building_mine <- ...

# try these operations
bld_itst <- st_intersection(building_mine, building_osm)
dif_osm_bld <- st_difference(building_osm, building_mine)
dif_bld_osm <- st_difference(building_mine, building_osm)
bld_union <- st_union(building_mine, building_osm)
bld_bbox <- st_bbox(bld_union) %>% 
  st_as_sfc() %>% st_sf()

# plot the results and learn each function
plot(bld_bbox$geometry, col = 'lightblue')
plot(building_osm$geometry, border = "black", add = T)
plot(building_mine$geometry, border = "brown", add = T)

plot(dif_osm_bld$geometry, col = 'red', add = T)
plot(dif_bld_osm$geometry, col = 'blue', add = T)
plot(bld_itst$geometry, col = 'yellow', add = T)

Homework

  • No homework
  • But please finish reading Unit2-Module1 and do all the practices (!!)