Punktinformationen aggregieren

Rafael Hologa
Dec 21, 2020 5 min read R
Image credit: Rafael Hologa

Bei der Verarbeitung von Geodaten werden oftmals Informationen, die auf der Individualebene in Form von Punktgeometrien vorliegen, für administrative Einheiten aggregiert. Eine solche Aggregation verursacht oft Artefakte, wie Zonierungs- und Skaleneffekte, sodass die Gefahr eines Ökologischen Fehlschlusses besteht. Um diese Effekte möglichst gering zuhalten, bieten sich Hexagone als Aggreagtionseinheit an.1

Problemstellung

Durch die Übertragung der Eigenschaften von Punktgeometrien auf künstliche Raumeinheiten - wie Gebietsflächen, regelmäßige Grids oder Hexagone - liegen die räumlichen Verteilungsmuster der Individualinformationen nicht mehr. Dabei entsteht ein neuer Datensatz, der suggeriert, dass innerhalb der Aggregationsflächen Eigenschaften ubiquitär und homogen vorliegen. Dies kann dazu führen, dass auf der Grundlage solcher Daten ein Ökologischer Fehlschluss gezogen wird.

Lösungsansätze

Damit bei der Aggregations von räumlichen Eigenschaften nur geringe Zonierungseffekte zum Tragen kommen, empfiehlt es sich keine Geometrien zu verwenden, die unterschiedlich große Flächen und verschiedene Formen aufweisen. Negative Skaleneffekte können abgeschwächt werden indem möglichst kleine Einheiten zur Aggregation verwendet werden, denn je kleiner die Flächen sind, desto positionsgetreuer bleiben auch die aggregierten Eigenschaften. Wegen ihrer großen Kompaktheit und da sie nebeneinander angeordnet eine Fläche lückenlos abdecken (ungleich ein Kreisen), bieten sich gleichmäßig angeordnete Hexagone als räumliche Aggregationseinheit an.

Punkte innerhalb Hexagone mit R aggregieren

Nötige R-Bbibliotheken laden

Zur Umsetzung des nachfolgenden Beispiels sind folgende R-Pakete nötig:

if (!require(sf)){install.packages("sf"); library(sf)}
if (!require(lwgeom)){install.packages("lwgeom"); library(lwgeom)}
if (!require(dplyr)){install.packages("dplyr"); library(dplyr)}
if (!require(ggplot2)){install.packages("ggplot2"); library(ggplot2)}
if (!require(units)) {install.packages("units"); library(units)}
if (!require(osmdata)){install.packages("osmdata"); library(osmdata)}

Beispieldaten laden

Die Individualinformationen werden beispielhaft von Punkten (‘points’) repräsentiert, die bestimmte Gebäudestandorte im Satdtgebiet von Freiburg im Breisgau (‘admins’) referenzieren.

# Reading layer 'statistische bezirke' as sf object
download.file("https://github.com/RafHo/teaching/raw/master/data/statistische_bezirke.zip", destfile = "statistische_bezirke.zip")
# unzip files
unzip("statistische_bezirke.zip")
# read shapefile
admins <- st_read("statistische_bezirke.shp")
## Reading layer `statistische_bezirke' from data source 
##   `V:\github\nils_website\blogdown\content\post\density\statistische_bezirke.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 169 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 400118.6 ymin: 5306198 xmax: 420233.5 ymax: 5324876
## Projected CRS: ETRS89 / UTM zone 32N
# retrieve points from osm-api
points <- opq(bbox = getbb("Freiburg, Germany")) %>% 
  add_osm_feature(key = "building", value= c("house", "residential",
                                             "residential", "apartments"))
points <- osmdata_sf(points)
points <- points$osm_points
#drop attributes
points <- points[,1]

### ALTERNATIVE APPROACH: random sample 10 000 points within the boundaries of 'admins'
#set.seed(1234)
#points <- st_sample(admins, 25, type = "random")

# transform points
points <- st_transform(points, crs=st_crs(admins))
# find points within polygons
points <- st_intersection(points, admins)
#drop attributes again
points <- points[,1]
# export points
st_write(points, "points.shp", update = T)
## Updating layer `points' to data source `points.shp' using driver `ESRI Shapefile'
## Updating existing layer points
## Writing 98649 features with 1 fields and geometry type Point.
# convert points to dataframe
points_df  <- cbind(points, st_coordinates(points))

#plot points
points_plot <- ggplot(data = admins )+
  geom_sf(fill = "white") +
  geom_point(data = points_df, aes(x = X, y = Y), size = 0.01) +
  xlab("") +
  ylab("") +
  theme_minimal()
points_plot

Mit Blick auf das Veteriliungsmuster der Punkte wird ersichtlich, dass sich diese Individualinformation sehr ungleich auf die unterschiedlichen statistische Bezirke verteilen.

Informationen aggregieren

Bei einer Aggregation der Punkte (Individualinformation) mittels der Polygonflächen (statistische Bezirke) wird ersichtlich, wieviele Punkte innerhalb der Polygone pro km2 auftreten. Hierfür werden die Punkte pro Polygon ausgezählt (topologischer Operator: Intersection) und auf die Gesamtfläche des jeweiligen Polygons normalisiert, da die Flächen sehr unterschiedlich groß sind.

# number of points within admin units
admins$number <- lengths(st_intersects(admins, points))

#define unit to km^2
admins$area <- set_units(st_area(admins), km^2)
# count density 
admins$count <- as.double(admins$number/admins$area)
# replace 0 values to NA
admins$count[admins$count == 0] <- NA

#find classes
brks <- seq(1, 16001, by=2000)

admins$size <- cut(admins$count,
                     breaks=brks,
                     labels=1:8,
                     include.lowest=TRUE)

admins$color <- cut(admins$count,
                     breaks=brks,
                     include.lowest=TRUE,
                     labels=c('#ffa83e','#ff8b24','#fd6a0b',
                              '#f04f09', '#d8382e', '#c62833','#AF1C43', '#8A1739'))

lab <- round(brks, digits = 0)
admins_ranges <- paste(lab,"-", lab[2:9])
admins_ranges <- admins_ranges[1:8]

admins$ranges <- cut(admins$count,
                     breaks=brks,
                     include.lowest=TRUE,
                     labels=admins_ranges)

#plot admin map
admins_plot <- ggplot()+
  geom_sf(data = admins, aes(fill = ranges), size = 0.1, colour = "white") +
  scale_fill_manual(values = levels(admins$color), name = expression("Punkte pro "~m^2)) +
  theme_minimal() 
admins_plot

Im Vergleich zu einer Aggregation der Punkte mittels gleichförmiger Hexagone wird ersichtlich, dass die Distanzen der Polygongrenzen zu den ursprünglichen Positionen der Individualdaten im Mittel geringer sind, als bei den ungleichförmigen und unterschiedlich großen Verwaltungsgrenzen. Je nach gewählter Skala, d.h Größe der Hexagone wirkt sich dieser Effekt stärker bzw. geringer aus.

#unit: meter
hexagons <- st_make_grid(admins, cellsize = 500, square=FALSE)
hexagons <- st_as_sf(hexagons)

# count points within admin units
hexagons$number <- lengths(st_intersects(hexagons, points))

#define unit to km^2
hexagons$area <- set_units(st_area(hexagons), km^2)
# count density 
hexagons$count <- as.double(hexagons$number/hexagons$area)
# replace 0 values to NA
hexagons$count[hexagons$count == 0] <- NA

#define classes
brks <- seq(1, 16001, by=2000)
hexagons$size <- cut(hexagons$count,
                     breaks=brks,
                     labels=1:8,
                     include.lowest=TRUE)

hexagons$color <- cut(hexagons$count,
                     breaks=brks,
                     include.lowest=TRUE,
                     labels=c('#ffa83e','#ff8b24','#fd6a0b',
                              '#f04f09', '#d8382e', '#c62833','#AF1C43', '#8A1739'))

lab <- round(brks, digits = 0)
hexagons_ranges <- paste(lab,"-", lab[2:9])
hexagons_ranges <- hexagons_ranges[1:8]
                                  

hexagons$ranges <- cut(hexagons$count,
                     breaks=brks,
                     include.lowest=TRUE,
                     labels=hexagons_ranges)

#plot grids/ hexagons map
hexagons_plot <- ggplot() + 
  geom_sf(data = hexagons, aes(fill = ranges), size = 0.1, colour = "white") +
  scale_fill_manual(values = levels(hexagons$color),  name = expression("Punkte pro "~m^2)) +
  theme_minimal() 
hexagons_plot

Dichte mittels Punktgrößen visualisieren

Um die Dichte der Punkte unabhängig von den zugrundeliegenden Flächen zu symbolisieren, kann der Ansatz auch dazu genutzt werden die Dichte anhand unterschiedlich großer Punkte zu visualisieren. Dafür werden die Zentroide der Hexgone ermittelt und die aggregierte Information zur Visualisierung der Punktgröße mitgeführt.

#create centroids
hexagons_cent <- st_centroid(hexagons)
hexagons_cent <- st_transform(hexagons_cent, crs=4326)
hexagons_cent <- hexagons_cent %>% filter(!count==0)
hexagons_cent$size <- as.numeric(hexagons_cent$size)

# convert points to dataframe
hexagons_cent_df  <- cbind(hexagons_cent, st_coordinates(hexagons_cent))
hexagons_cent_df$size <- as.numeric(hexagons_cent_df$size/1.5)

plot_hex_df <- ggplot() +
  geom_point(data = hexagons_cent_df, aes(color = ranges, x=X, y=Y), size = hexagons_cent_df$size) +
  scale_color_manual(values = levels(hexagons_cent_df$color),
                     name = expression("Punkte pro "~m^2)) +
  coord_quickmap() + #to scale the x and y axis equally by long and lat values
  theme_minimal() +
  xlab("Longitude") +
  ylab("Latitude") #+ theme(legend.position = "none")
plot_hex_df

  1. S. u.a.: Birch, C.P.D.; Oom, S.P.; Beecham, J.A. Rectangular and hexagonal grids used for observation, experiment and simulation in ecology. Ecol. Modell. 2007, 206, 347–359, doi:10.1016/j.ecolmodel.2007.03.041.↩︎

Hot-Spots Punkte GIS Aggregation Dichte Zonierung Skala
Avatar
Rafael Hologa
Research associate

My research interests are the Crowdsourcing of Volunteered Geographic Informations (VGI) & Harvesting Ambient Geospaial Informations (AGI) relating to perceptions of weather parameters in the context of open science and citizen science initiatives.

Related