Skip to content

Commit

Permalink
ex1
Browse files Browse the repository at this point in the history
  • Loading branch information
zumbov2 committed Mar 19, 2021
1 parent ac30295 commit c905bfd
Show file tree
Hide file tree
Showing 4 changed files with 234 additions and 1 deletion.
9 changes: 8 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
[![Build Status](https://travis-ci.org/zumbov2/swissgd.svg?branch=master)](https://travis-ci.org/zumbov2/swissgd)

# swissgd
This R package is an interface to parts of the [Geo-Information Platform of the Swiss Confederation](https://www.geo.admin.ch/en/geo-services/geo-services/download-services.html) and **work-in-progress**. It provides functions to search and download data from [data.geo.admin.ch](https://data.geo.admin.ch/) and wraps around the [Spatial Temporal Asset Catalog (STAC) API](https://data.geo.admin.ch/api/stac/v0.9/).
This R package is an interface to parts of the [Geo-Information Platform of the Swiss Confederation](https://www.geo.admin.ch/en/geo-services/geo-services/download-services.html) and **work-in-progress**. It provides functions to search and download geodata from [data.geo.admin.ch](https://data.geo.admin.ch/) and wraps around the [Spatial Temporal Asset Catalog (STAC) API](https://data.geo.admin.ch/api/stac/v0.9/).

The acquisition and use of data or services is free of charge, subject to the provisions on fair use. For more information, please see the [Terms of Use](https://www.geo.admin.ch/en/geo-services/geo-services/terms-of-use.html).

Expand Down Expand Up @@ -132,7 +132,14 @@ res <- swissgd::get_stac_assets(
swissgd::download_stac_assets(res)
```

# Examples

## Spatial distribution of place name suffixes
*Datasets:* ch.swisstopo.swissnames3d, ch.swisstopo.swissboundaries3d-land-flaeche.fill
*Packages:* `sf`, `raster`, `btb`, `ggplot2` etc.

<img src="https://github.com/zumbov2/swissgd/blob/master/examples/ex1_ikon.png" width="300">
<img src="https://github.com/zumbov2/swissgd/blob/master/examples/ex1_berg_mont.png" width="300">



Binary file added examples/ex1_berg_mont.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added examples/ex1_ikon.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
226 changes: 226 additions & 0 deletions examples/ex1_swissnames.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,226 @@
# Load Packages --------------------------------------------------------------
install.packages("devtools")
devtools::install_github("zumbov2/swissgd")

library(swissgd)
library(sf)
library(raster)
library(dplyr)
library(stringr)
library(purrr)
library(fasterize)
library(btb)
library(tidyr)
library(ggplot2)
library(hrbrthemes)

# Load Data ------------------------------------------------------------------

# Search dataset with place names
swissgd::search_geodata("Names")

# Download
download_geodata("ch.swisstopo.swissnames3d$")

# Unpack the zip folder inside the downloaded folder
unzip("ch.swisstopo.swissnames3d/swissNAMES3D_LV03.zip")

# Load the shapefile dataset and filter for places
ply <- sf::st_read("swissNAMES3D_LV03/shp_LV03_LN02/swissNAMES3D_PLY.shp")

places <- ply %>%
filter(OBJEKTART == "Ort") %>%
filter(STATUS == "offiziell")

# Place name endings ---------------------------------------------------------

# Get 2 to 7 letter place name endings
endings <- places %>%
mutate(chars = nchar(NAME)) %>%
mutate(
ending2 = stringr::str_sub(NAME, chars-1, chars),
ending3 = stringr::str_sub(NAME, chars-2, chars),
ending4 = stringr::str_sub(NAME, chars-3, chars),
ending5 = stringr::str_sub(NAME, chars-4, chars),
ending6 = stringr::str_sub(NAME, chars-5, chars),
ending7 = stringr::str_sub(NAME, chars-6, chars)
) %>%
select(starts_with("ending")) %>%
pivot_longer(starts_with("ending")) %>%
group_by(name, value) %>%
count() %>%
ungroup() %>%
arrange(desc(n))

# Select the most common ones manually
endings_selection <- c(
"hof", "berg", "matt", "egg", "bach", "wil", "weid", "haus", "au", "moos",
"holz", "acker", "acher", "boden", "tal", "bühl", "rüti", "hus", "hüsli", "büel",
"feld", "wald", "schwand", "rain", "ried", "stein", "ikon", "halde", "loch", "graben",
"halden", "mühle", "burg", "alp", "mont", "hubel"
)

# Detect endings in place names...
ending_detect <- purrr::map_dfc(
paste0(endings_selection, "$"),
function(x) str_detect(tolower(orte$NAME), x)
)

# ...bind to dataset and reduce geoms to points
names(ending_detect) <- endings_selection
places <- bind_cols(places, ending_detect) %>%
st_point_on_surface()

# Get national border as bounding --------------------------------------------

# Search
swissgd::search_geodata("boundaries")

# Download
swissgd::download_geodata("ch.swisstopo.swissboundaries3d-land-flaeche.fill")

# Readme says we find the data here
folder <- "swissboundaries3d.zip"
download.file(
"https://data.geo.admin.ch/ch.swisstopo.swissboundaries3d-gemeinde-flaeche.fill/shp/2056/ch.swisstopo.swissboundaries3d-gemeinde-flaeche.fill.zip",
folder
)
unzip(folder, exdir = "swissboundaries3d")
file.remove(folder)

# Load data
ch <- sf::st_read("swissboundaries3d/swissBOUNDARIES3D_1_3_TLM_LANDESGEBIET.shp")
ch <- ch %>% filter(NAME == "Schweiz")

# Transforming points pattern to continuous coverage -------------------------------------

# Strongly inspired by
# https://www.r-bloggers.com/2019/05/kernel-spatial-smoothing-transforming-points-pattern-to-continuous-coverage/
get_spatial_smooting <- function(ending, x, bounding, res, bw, crs = 2056) {

cat(ending, "\n")

# Bounding
bx <- sf::st_bbox(bounding)

# Make grid
bounding_coords <- bounding %>%
sf::st_make_grid(
cellsize = res,
offset = c(
plyr::round_any(bx[1] - bw, res, floor),
plyr::round_any(bx[2] - bw, res, floor)),
what = "centers"
) %>%
sf::st_sf() %>%
sf::st_join(bounding, join = st_intersects, left = F) %>%
sf::st_coordinates() %>%
tibble::as_tibble() %>%
dplyr::select(x = X, y = Y)

# Compute Kernel
kernel <- x %>%
cbind(., st_coordinates(.)) %>%
sf::st_set_geometry(NULL) %>%
dplyr::select(x = X, y = Y, ending) %>%
btb::kernelSmoothing(
dfObservations = .,
sEPSG = crs,
iCellSize = res,
iBandwidth = bw,
vQuantiles = NULL,
dfCentroids = bounding_coords
)

names(kernel)[3] <- "density"
kernel$ending <- paste0("-", ending)

return(kernel)

}

# Spatial Kernel Smooting for all endings
endings_kernels <- purrr::map_dfr(
endings_selection,
get_spatial_smooting,
x = places,
bounding = ch,
res = 1000,
bw = 20000
)

# Visualize ----------------------------------------------------------------

# ikon
p1 <- endings_kernels %>%
# filter(ending %in% c("-hof", "-matt", "-wil", "-ikon", "-berg", "-mont")) %>%
# mutate(ending = factor(ending, levels = c("-hof", "-matt", "-wil", "-ikon", "-berg", "-mont"))) %>%
filter(ending == "-ikon") %>%
# mutate(ending = factor(ending, levels = c("-hof", "-matt", "-wil", "-ikon", "-berg", "-mont"))) %>%
group_by(ending) %>%
mutate(density = density / max(density)) %>%
ungroup() %>%
ggplot() +
geom_sf(fill = "grey90", color = NA) +
geom_sf(aes(fill = density, alpha = density), color = NA) +
# facet_wrap(.~ending) +
scale_fill_viridis_c(option = "magma", direction = -1, na.value = "white") +
scale_alpha_continuous(range = c(0, 1)) +
labs(
title = "Distribution of the place name suffix -ikon",
subtitle = "Spatial kernel density estimation",
caption = "Datasets: ch.swisstopo.swissnames3d, ch.swisstopo.swissboundaries3d
Method: btb::kernelSmoothing with a grid resolution of 1km and a bandwidth of 20km"
) +
theme_ipsum_rc() +
theme(
legend.position = "none",
plot.background = element_rect(fill = "white", color = NA),
legend.background = element_rect(fill = "white", color = NA),
panel.background = element_rect(fill = "white", color = NA),
axis.text = element_text(color = "white"),
axis.ticks = element_line(color = "white"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.caption = element_text(hjust = 0)
)

ggsave("ikon.png", p1, width = 11.1/1.5, height = 8.33/1.5, dpi = 500)

# Suffix for mountain in german and french
p2 <- endings_kernels %>%
filter(ending %in% c("-berg", "-mont")) %>%
rename("Suffix" = "ending") %>%
group_by(Suffix) %>%
mutate(density = density / max(density)) %>%
ungroup() %>%
ggplot() +
geom_sf(fill = "grey90", color = NA) +
geom_sf(aes(fill = Suffix, alpha = density), color = NA) +
scale_alpha_continuous(range = c(0, 1), guide = FALSE) +
scale_fill_manual(values = c("#377eb8", "#4daf4a")) +
labs(
title = "Distribution of the place name suffixes -berg and -mont",
subtitle = "Spatial kernel density estimation",
caption = "Datasets: ch.swisstopo.swissnames3d, ch.swisstopo.swissboundaries3d
Method: btb::kernelSmoothing with a grid resolution of 1km and a bandwidth of 20km"
) +
theme_ipsum_rc() +
theme(
legend.position = "right",
legend.justification = c(0, 1),
plot.background = element_rect(fill = "white", color = NA),
legend.background = element_rect(fill = "white", color = NA),
panel.background = element_rect(fill = "white", color = NA),
axis.text = element_text(color = "white"),
axis.ticks = element_line(color = "white"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.caption = element_text(hjust = 0)
)

ggsave("berg_v_mont.png", p2, width = 11.1/1.5, height = 8.33/1.5, dpi = 500)

0 comments on commit c905bfd

Please sign in to comment.