Aggregation Demo

Lade Postleitzahl-Gebiete und R-Erweiterungen

library(rgdal)

## Loading required package: sp

## rgdal: version: 1.3-1, (SVN revision 747)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: /usr/share/gdal/2.2
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.2-7

library(lattice)
library(latticeExtra)

## Loading required package: RColorBrewer

library(rgeos)

## rgeos version: 0.3-27, (SVN revision 560)
##  GEOS runtime version: 3.6.2-CAPI-1.10.2 4d2925d6 
##  Linking to sp version: 1.2-7 
##  Polygon checking: TRUE

plz_poly <- readOGR("postleitzahlengebiete_mecklenburg-vorpommern.kml")

## OGR data source with driver: LIBKML 
## Source: "/home/matthias/plz_aggregation/postleitzahlengebiete_mecklenburg-vorpommern.kml", layer: "postleitzahlengebiete_mecklenburg-vorpommern"
## with 190 features
## It has 13 fields

plz <- as.numeric(as.character(plz_poly$postleitzahl))
# wähle alle PLZ in Rostock
rostock_poly <- plz_poly[plz %in% c(18055, 18057, 18059, 18069, 18106, 18107, 18109, 18119, 18146, 18147, 18181, 18182),]
## Projiziere Koordinaten nach Gauß-Krüger (Zone 3)
proj4string(plz_poly)

## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

rostock_poly_gk3 <- spTransform(rostock_poly, CRSobj = CRS("+init=epsg:31467"))

Erzeuge Zufallsdaten für 2x2 km Raster

# Ereuge Raster
rostock_grid <- makegrid(rostock_poly_gk3, cellsize = 2000)
gridded(rostock_grid) <- ~x1+x2
fullgrid(rostock_grid)<- TRUE
proj4string(rostock_grid)<- proj4string(rostock_poly_gk3)
#Erzeuge Attribut-Daten
dummydata <- data.frame(values = runif(length(rostock_grid)))
rostock_grid <- SpatialGridDataFrame(rostock_grid, dummydata)

Erzeuge aggregierte Daten

aggregated_mean <- over(rostock_poly_gk3, rostock_grid, fn = mean)
aggregated_sd <- over(rostock_poly_gk3, rostock_grid, fn = sd)
result_poly <- SpatialPolygonsDataFrame(rostock_poly_gk3, data.frame(mean = aggregated_mean, sd = aggregated_sd))

Ergebnis

rostock_grid_poly <- as(rostock_grid, "SpatialPolygons")
p <- spplot(result_poly, names.attr=c("Mittelwert","Standardabweichung"))
p +
  layer(sp.polygons(rostock_grid_poly)) +
  layer(sp.text(coordinates(rostock_grid), txt=round(rostock_grid$values,digits = 2), cex=0.5))