Week 12: Working with OpenStreetMap Data and Data in Raster Format

PPOL 6805: GIS for Spatial Data Science

Christy Hsu

th1010@georgetown.edu

2026-11-13

R Coding Workshop for GIS: 10th Meeting

Outline

• Open Street Map
• Project checklist
• Raster Data

Open Street Map

1. use place name to search for bounding box

?getbb

taipei_bb_results <- getbb(place_name = 'taipei', format_out='sf_polygon')
taipei_bb_results

place_id	licence	osm_type	osm_id	lat	lon	class	type	place_rank	importance	addresstype	name	display_name	geometry
231159944	Data © OpenStreetMap contributors, ODbL 1.0. http://osm.org/copyright	relation	1293250	25.0375198	121.5636796	boundary	administrative	7	0.7098064	city	臺北市	臺北市, 臺灣	POLYGON ((121.4571 25.10798…
416923803	Data © OpenStreetMap contributors, ODbL 1.0. http://osm.org/copyright	relation	1527220	25.0119970	121.4656619	boundary	administrative	8	0.6239171	city	新北市	新北市, 臺灣	MULTIPOLYGON (((121.2826 25…

# tpc_bb <- taipei_bb_results |> filter(place_id == 213710844)
tpc_bb <- taipei_bb_results |> filter(place_id == taipei_bb_results$place_id[1])
tpc_bb |> mapview()

# ntpc_bb <- taipei_bb_results |> filter(place_id == 212454050)
ntpc_bb <- taipei_bb_results |> filter(place_id == taipei_bb_results$place_id[2])
ntpc_bb |> mapview()

osmdata

2. create a query

# assign("has_internet_via_proxy", TRUE, environment(curl::has_internet))

tpc_mosque_results <- tpc_bb |>
  opq() |> 
  add_osm_feature(key='building', value = 'mosque') |>
  osmdata_sf()
tpc_mosque_results

# saveRDS(
#   tpc_mosque_results,
#   'data/tpc-mosque-results.rds'
# )
tpc_mosque_results <- readRDS('data/tpc-mosque-results.rds')

3. leaflet visualization

leaflet() |>
  addProviderTiles('CartoDB.Positron') |>
  addPolygons(
    data = tpc_mosque_results$osm_polygons,
    color = "purple"
  )

tpc_mosque_polygon_sf <- tpc_mosque_results$osm_polygons

osmdata

kentucky_bb <- getbb('kentucky, USA')
kentucky_bb

        min       max
x -89.57151 -81.96454
y  36.49712  39.14780

1. Bounding Box

ky_results <- kentucky_bb |>
  opq() |>
  add_osm_feature(key = 'boundary', value = 'administrative') |>
  add_osm_feature(key = 'admin_level', value = '4') |>
  add_osm_feature(key = 'name', value = 'Kentucky') |>
  osmdata_sf()
ky_results

# ky_results |> saveRDS('data/ky-result.rds')

ky_results <- readRDS('data/ky-result.rds')
ky_sf <- ky_results$osm_multipolygons
ky_vars <- c("osm_id", "name", "ISO3166-2", "admin_level", "border_type", 
"boundary", "nickname", "official_name", "ref", "ref:USCG", "ref:USPS", "ref:fips", 
"short_name", "source:short_name", "type", 
"wikidata", "wikipedia", "geometry")

ky_sf <- ky_sf |> dplyr::select(all_of(ky_vars))
ky_sf |> mapview()

ky_bookcase_results <- kentucky_bb |>
  opq() |>
  add_osm_feature(key='amenity', value = 'public_bookcase') |>
  osmdata_sf()
ky_bookcase_results

# ky_bookcase_results |> saveRDS('data/ky-bookcase-results.rds')

ky_bookcase_results <- readRDS('data/ky-bookcase-results.rds')

ky_bookcase_sf <- ky_bookcase_results$osm_points

# ky_convex_sf <- ky_sf |> st_convex_hull() |> st_buffer()
# ky_convex_sf |> mapview()
ky_bbox_sf <- ky_sf |> st_bbox() |> st_as_sfc() |> st_as_sf(crs = 4326)
ky_bbox_sf |> mapview()

Tessellation

library(tigris)

uac_sf <- urban_areas(year = 2020, progress_bar = FALSE)

uac_sf |> st_crs() |> print()

Coordinate Reference System:
  User input: NAD83 
  wkt:
GEOGCRS["NAD83",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4269]]

ky_sf |> st_crs() |> print()

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

uac_sf <- uac_sf |> st_transform(4326)

ky_uac_sf <- uac_sf[ky_sf,]
ky_urban_bookcase_sf <- ky_bookcase_sf[ky_uac_sf,]
ky_urban_bookcase_sf |> dim()

[1] 68 24

mapview(ky_sf, col.region = 'lightblue', alpha.regions = 0.6, map.types = 'CartoDB.Positron') +
  mapview(ky_uac_sf, col.regions = 'yellow') + 
  mapview(ky_urban_bookcase_sf, cex = 4)

Observation Window

kybb_uac_sf <- uac_sf[ky_bbox_sf,]
kybb_uac_sf |> mapview()

ky_bookcase_sf[kybb_uac_sf,] |> dim()

[1] 75 24

ppp

ky_owin <- ky_bbox_sf |> st_transform(3857) |> as.owin()

ky_urban_owin <- kybb_uac_sf |>
  st_simplify(dTolerance = 200) |>
  st_intersection(ky_bbox_sf) |>
  st_transform(3857) |>
  as.owin()
# plot(ky_urban_owin)

ky_bookcase_ppp <- ky_bookcase_sf |>
  st_as_sfc() |>
  st_transform(3857) |>
  as.ppp(W = ky_owin)
plot(ky_bookcase_ppp)

Dividing the observation window

urban_rural_tess <- tess(
  tiles = list(
    urban = ky_urban_owin,
    rural = setminus.owin(ky_owin, ky_urban_owin)
  )
)

urban_rural_tess |> class()

[1] "tess" "list"

Simulation

set.seed(6805)
gen_sim_ppp <- function() {
  prot_sim <- spatstat.random::rpoint(
    n = nrow(ky_bookcase_sf),
    w = ky_owin
  )
  return(prot_sim)
}
sim_prot_ppp <- gen_sim_ppp()
plot(sim_prot_ppp)

compute_tile_counts <- function(ppp_obj) {
  class_list <- tileindex(
    ppp_obj,
    Z = urban_rural_tess
  )
  count_tb <- table(class_list)
  count_df <- tibble(
    urban = count_tb[1],
    rural = count_tb[2]
  )

  return(count_df)
}

compute_tile_counts(sim_prot_ppp)

urban	rural
3	82

compute_tile_counts(ky_bookcase_ppp)

urban	rural
73	12

Simulation

set.seed(6805)
gen_sims_ppp <- function(num_sims = 999) {
  bkcase_sims <- spatstat.random::rpoint(
    n = nrow(ky_bookcase_sf),
    w = ky_owin,
    nsim = num_sims
  )
  return(bkcase_sims)
}
bkcase_sims_list <- gen_sims_ppp()

# compute_tile_counts <- function(ppp_obj) {
#   class_list <- tileindex(
#     ppp_obj,
#     Z = urban_rural_tess
#   )
#   counts <- class_list |>
#     table() |>
#     as.vector()
#   names(counts) <- c('urban', 'rural')

#   return(counts)
# }

sims_tess_counts <- lapply(
  X = bkcase_sims_list,
  FUN = compute_tile_counts
)
sim_df <- sims_tess_counts |> bind_rows()
sim_df |> slice_sample(n = 4)

urban	rural
4	81
2	83
3	82
1	84

obs_df <- compute_tile_counts(ky_bookcase_ppp)
obs_df

urban	rural
73	12

bind_df <- bind_rows(sim_df, obs_df)
bind_df |> dim()

[1] 1000    2

bind_df |>
  ggplot(aes(x = rural)) +
  geom_density(fill = palette4[2], alpha = 0.7) +
  geom_vline(
    xintercept = 9,
    linetype = 'dashed',
    color = palette4[5]
  ) +
    theme_classic()

Raster

pop_raster <- raster(
  'data/gpw_v4_population_count_adjusted_to_2015_unwpp_country_totals_rev11_2020_2pt5_min.tif'
)
pop_raster

class      : RasterLayer 
dimensions : 4320, 8640, 37324800  (nrow, ncol, ncell)
resolution : 0.04166667, 0.04166667  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : gpw_v4_population_count_adjusted_to_2015_unwpp_country_totals_rev11_2020_2pt5_min.tif 
names      : gpw_v4_population_count_adjusted_to_2015_unwpp_country_totals_rev11_2020_2pt5_min 

bbox_extent <- extent(-89.57151, -81.96454, 36.49712, 39.14780)
pop_raster_sub <- crop(pop_raster, bbox_extent)

pop_poly <- pop_raster_sub |> raster::rasterToPolygons(n = 8, na.rm = TRUE)
pop_poly

class       : SpatialPolygonsDataFrame 
features    : 11712 
extent      : -89.58333, -81.95833, 36.5, 39.16667  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       : gpw_v4_population_count_adjusted_to_2015_unwpp_country_totals_rev11_2020_2pt5_min 
min values  :                                                                                 0 
max values  :                                                                     37111.3984375 

pop_sf <- st_as_sf(pop_poly)
pop_sf |> glimpse()

Rows: 11,712
Columns: 2
$ gpw_v4_population_count_adjusted_to_2015_unwpp_country_totals_rev11_2020_2pt5_min <dbl> …
$ geometry                                                                          <POLYGON [°]> …

pop_sf <- pop_sf |> rename(pop_count = gpw_v4_population_count_adjusted_to_2015_unwpp_country_totals_rev11_2020_2pt5_min)

gpw_map <- pop_sf |>
  ggplot(aes(fill = pop_count)) +
  geom_sf(color = scales::alpha("white",0)) +
  scale_fill_viridis_c(trans = 'pseudo_log') + 
  theme_minimal() +
  labs(
    title = "Population Counts",
    subtitle = 'GPW V.4 2020-07-01',
    fill = "Population Count"
    ) + 
  theme(
    plot.title = element_text(hjust=0.5),
    plot.subtitle = element_text(hjust=0.5)
    )

ggsave('image/gpw-raster.png', plot = gpw_map)

gpw_map