Worldwide Berlin is famously known for its nightlife. The data set we will work with stores the geographical location of clubs and bars located in Berlin. The data source is gathered from OpenStreetMap (OSM). The data was downloaded from GEOFABRIK on August 27, 2022, and contains OpenStreetMap data as of August 26, 2022 (see here).
We download the data and read the
gis_osm_pois_free_1.shp file, which corresponds to
points of interests in Berlin, using the st_read()
function from the sf package.
## build a temporary folder on disk
temp <- tempfile()
download_url <- "https://userpage.fu-berlin.de/soga/data/raw-data/spatial/"
zipfile <- "berlin-latest-free.shp.zip"
## download the file
download.file(paste0(download_url, zipfile), temp, mode = "wb")
## unzip the file(s)
unzip(temp)
## close file connection
unlink(temp)library(rgeos)
library(sf)
library(rgdal)
## read in the data
berlin_features <- st_read("./gis_osm_pois_free_1.shp", quiet = TRUE)Applying the str() function gives us more information
about the structure of the data set.
str(berlin_features)## Classes 'sf' and 'data.frame': 83123 obs. of 5 variables:
## $ osm_id : chr "16541597" "26735749" "26735753" "26735759" ...
## $ code : int 2907 2301 2006 2301 2301 2701 2307 2031 2724 2907 ...
## $ fclass : chr "camera_surveillance" "restaurant" "telephone" "restaurant" ...
## $ name : chr "Aral" "Aida" NA "Madame Ngo" ...
## $ geometry:sfc_POINT of length 83123; first list element: 'XY' num 13.3 52.5
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA
## ..- attr(*, "names")= chr [1:4] "osm_id" "code" "fclass" "name"The data set of class sf and data.frame
contains 83123 records (denoted as observations) represented as
rows. The 5 variables are represented by columns, with the first one
containing the osm_id number. From the fifth variable
geometry we know that the data set stores information about
83123 two-dimensional ('XY') surface points. The variable
fclass contains the category related to a point, such as
“restaurant” or “camera_surveillance”.
By applying the unique() function we get an overview of
all the different categories represented in the data set.
unique(berlin_features$fclass)## [1] "camera_surveillance" "restaurant" "telephone"
## [4] "tourist_info" "biergarten" "recycling_glass"
## [7] "memorial" "fire_station" "artwork"
## [10] "bank" "post_office" "toilet"
## [13] "cinema" "library" "police"
## [16] "post_box" "chemist" "motel"
## [19] "kiosk" "arts_centre" "cafe"
## [22] "fast_food" "car_sharing" "viewpoint"
## [25] "pub" "hotel" "supermarket"
## [28] "bakery" "pharmacy" "atm"
## [31] "convenience" "theatre" "hostel"
## [34] "car_rental" "gift_shop" "museum"
## [37] "guesthouse" "shelter" "stationery"
## [40] "bar" "travel_agent" "embassy"
## [43] "butcher" "beverages" "pitch"
## [46] "outdoor_shop" "playground" "kindergarten"
## [49] "attraction" "school" "monument"
## [52] "hunting_stand" "sports_centre" "tower"
## [55] "bench" "recycling_clothes" "optician"
## [58] "fountain" "nightclub" "drinking_water"
## [61] "computer_shop" "furniture_shop" "recycling"
## [64] "community_centre" "bicycle_shop" "doctors"
## [67] "veterinary" "bookshop" "clothes"
## [70] "florist" "doityourself" "laundry"
## [73] "comms_tower" "hairdresser" "university"
## [76] "video_shop" "picnic_site" "car_dealership"
## [79] "beauty_shop" "dentist" "vending_any"
## [82] "vending_machine" "mobile_phone_shop" "track"
## [85] "newsagent" "shoe_shop" "waste_basket"
## [88] "toy_shop" "town_hall" "sports_shop"
## [91] "vending_cigarette" "greengrocer" "water_well"
## [94] "mall" "jeweller" "bicycle_rental"
## [97] "department_store" "garden_centre" "water_tower"
## [100] "vending_parking" "theme_park" "nursing_home"
## [103] "car_wash" "graveyard" "college"
## [106] "recycling_paper" "water_works" "chalet"
## [109] "ruins" "clinic" "zoo"
## [112] "market_place" "general" "battlefield"
## [115] "stadium" "archaeological" "hospital"
## [118] "camp_site" "swimming_pool" "prison"
## [121] "observation_tower" "wayside_shrine" "water_mill"
## [124] "caravan_site" "golf_course" "courthouse"
## [127] "food_court" "dog_park" "wayside_cross"
## [130] "park"There are 130 different categories.
Since we want to have a closer look at Berlin’s nightlife, we subset
our data to include only the category nightclub.
berlin_locations <-
berlin_features[
berlin_features$fclass == "nightclub",
c("fclass", "name")
]Now, we also want to include bars in our data set.
Exercise: Subset the data set
berlin_featuressuch that it only includes the categoriesnightclubandbarby using the boolean operatorOR.
## Your code here...
berlin_locations <- NULLShow code
berlin_locations <-
berlin_features[
berlin_features$fclass == "nightclub" | berlin_features$fclass == "bar", # "OR" and "|" are equivalent
c("fclass", "name")
]
We rename the column fclass into loc and
plot the relative frequency of the categories in our data set by
combining the table(), the prop.table() and
the barplot() function.
# rename column from fclass to loc
colnames(berlin_locations)[1] <- "loc"
# plot
barplot(prop.table(table(berlin_locations$loc)))
As expected, there are significantly more locations denoted as bar in our data set than locations denoted as nightclub. Let us compare the absolute figures.
Exercise: Determine the absolute numbers of the locations
barandnightclubof our data set.
## Your code here...Show code
table(berlin_locations$loc)##
## bar nightclub
## 807 131
We visualize the data with the mapview package.
library(mapview)
mapview(berlin_locations,
zcol = "loc",
label = berlin_locations$name,
legend = TRUE
)Finally, we project the berlin_locations data set to the
European Terrestrial Reference System 1989 (ETRS89/UTM
zone 32N).
Exercise: Transform the coordinate reference system of the data set
berlin_locationsinto ERTS89/UTM zone 32N (EPSG 25832). Please visit the section about coordinate reference systems and transformations on the previous page if you want to look it up again.
## Your code here...
berlin_locations <- NULLShow code
berlin_locations <- st_transform(berlin_locations, 25832)
st_geometry(berlin_locations)## Geometry set for 938 features
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 780782 ymin: 5813041 xmax: 814411.7 ymax: 5841533
## Projected CRS: ETRS89 / UTM zone 32N
## First 5 geometries:## POINT (793357 5832558)## POINT (797981.4 5827060)## POINT (797927.1 5827157)## POINT (792808.9 5832719)## POINT (799189.1 5830723)
Regular spaced and random data
In order to contrast particular statistical methods for spatial point
pattern analysis we construct two artificial data sets. One data set,
berlin_random, consists of point data randomly distributed
over the area of Berlin. The other data set, berlin_regular
consists of point data that is fairly regular distributed across the
area of Berlin.
Before we start with the data generation we have to set the
geographical extend. Hence, we load the the sf object
containing Berlin’s district borders (berlin_sf).
load(url("https://userpage.fu-berlin.de/soga/data/r-data/berlin_district.RData"))
plot(berlin_sf, main = "")
We take the bounding box specifications of the berlin_sf
variable to set the geographic extend for the data generation
process.
xmin <- as.vector(st_bbox(berlin_sf)[1])
ymin <- as.vector(st_bbox(berlin_sf)[2])
xmax <- as.vector(st_bbox(berlin_sf)[3])
ymax <- as.vector(st_bbox(berlin_sf)[4])
Random points
First, we generate random uniform numbers using the
runif() function.
set.seed(1111) # set seed for reproducibility
## generate random uniform numbers
x_random <- runif(
n = nrow(berlin_locations) * 0.4,
min = xmin,
max = xmax
)
y_random <- runif(
n = nrow(berlin_locations) * 0.4,
min = ymin,
max = ymax
)Then we combine the random number vector data in a
data.frame object and plot it together with the district
border of Berlin.
## add to data frame
berlin_random <- data.frame("x" = x_random, "y" = y_random)
plot(
x = berlin_random$x,
y = berlin_random$y,
main = "Random Points",
xlab = "x", ylab = "y", cex = .5,
asp = TRUE
)
plot(berlin_sf, add = TRUE, border = "blue", col = NA)
Regular points
First, we generate a sequence of numbers using the seq()
function.
x_regular <- seq(
from = xmin, to = xmax,
length.out = sqrt(nrow(berlin_locations)) * 0.35
)
y_regular <- seq(
from = ymin, to = ymax,
length = sqrt(nrow(berlin_locations)) * 0.35
)Then we create a data frame from all combinations of variables using
the expand.grid() function.
regular <- expand.grid(x_regular, y_regular)
plot(regular, asp = TRUE)
Finally, we add some noise to the points using the
jitter() function, combine the data into a
data.frame object, berlin_regular, and plot it
together with the district border of Berlin.
x_regular <- jitter(regular[, 1], factor = 1)
y_regular <- jitter(regular[, 2], factor = 1)
berlin_regular <- data.frame("x" = x_regular, "y" = y_regular)
plot(
x = berlin_regular$x,
y = berlin_regular$y,
main = "Regular Points",
xlab = "x", ylab = "y", cex = .5,
asp = TRUE
)
plot(berlin_sf, add = TRUE, border = "blue", col = NA)
Creating spatstat objects
In spatstat a spatial point pattern in two-dimensional
space is stored as an object of class ppp. Such a point
pattern object contains spatial coordinates of the
points, marks attached to them (if any), a
window in which the points were observed and optionally
the name of the unit of length for the spatial
coordinates.
A ppp object is created using the ppp()
function. The argument window must be provided a window
object of class owin. Including information about the units
of length in which the \(x\) and \(y\) coordinates are recorded, enables the
spatstat package to print better reports and to annotate
the axes in plots.
library(spatstat)
library(maptools)For the purpose of this tutorial we will build three ppp
objects, one for the berlin_locations data set which we
denote as ppp_locations, one for the
berlin_random data set (ppp_random) and one
for the berlin_regular data set
(ppp_regular).
All of the three data sets share the same unit of length
(meter) in which the \(x\) and
\(y\) coordinates are recorded as well
as the window in which the points were observed (the bounding box of the
Berlin district border berlin_sf).
Coordinate units
We start with assigning a placeholder variable which holds the string representation of the coordinate’s units.
coordinate_units <- "meters"
Window
An object of class owin (observation window)
represents a region or window in a two-dimensional space. The window may
be a rectangle, a polygon or polygons, with polygonal holes, or an
irregular shape represented by a binary pixel image mask.
Recall that berlin_sf is of class sf.
class(berlin_sf)## [1] "sf" "data.frame"To set the window to the polygon berlin_sf, we first
need to coerce berlin_sf into a Spatial
object, which then is coerced into an owin object. This may
sound complicated, but by using the as() function in
sequence the objective becomes fairly simple.
class(as(berlin_sf, "Spatial"))## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"class(as(as(berlin_sf, "Spatial"), "owin"))## [1] "owin"Here we go, done! Now we assign the owin object to a
variable denoted as win.
win <- as(as(berlin_sf, "Spatial"), "owin")
win## window: polygonal boundary
## enclosing rectangle: [777974.1, 823510.5] x [5808837, 5845580] units
The ppp object
Finally, we construct ppp objects using the
ppp() function from our three different data sets
berlin_locations, berlin_regular and
berlin_random. Recall that berlin_locations is
an sf object, whereas berlin_regular and
berlin_random are data.frame objects.
class(berlin_locations)## [1] "sf" "data.frame"class(berlin_regular)## [1] "data.frame"class(berlin_random)## [1] "data.frame"For the sf object, we use the
st_coordinates() function to extract the \(xy\)-coordinates in form of a
matrix object. Let us review the first five point
coordinates in berlin_locations in matrix format.
st_coordinates(berlin_locations)[1:5, ]## X Y
## [1,] 793357.0 5832558
## [2,] 797981.4 5827060
## [3,] 797927.1 5827157
## [4,] 792808.9 5832719
## [5,] 799189.1 5830723Now, let us construct the ppp object.
ppp_locations <- ppp(
x = st_coordinates(berlin_locations)[, 1],
y = st_coordinates(berlin_locations)[, 2],
window = win
)## Warning: 4 points were rejected as lying outside the specified windowNow, we assign the name of the unit of length for the spatial
coordinates to the ppp object and plot it.
unitname(ppp_locations) <- coordinate_units
plot(ppp_locations)## Warning in plot.ppp(ppp_locations): 4 illegal points also plotted
Note that we received a warning that 4 points were rejected as
lying outside the specified window. Do not worry, the rejected
points in a point pattern will be ignored by all other functions except
for plot. They are still accessible and stored as an
attribute of the point pattern called rejects.
It is not very easy to spot the rejected points in the plot above,
thus let us write some code to highlight the rejected points. We use the
gIntersects() function from the rgeos
package.
sp_locations <- as(berlin_locations, "Spatial")
sp_border <- as(berlin_sf, "Spatial")
i <- gIntersects(sp_locations, sp_border, byid = TRUE)## Warning in RGEOSBinPredFunc(spgeom1, spgeom2, byid, func): spgeom1 and spgeom2
## have different proj4 stringsplot(sp_border)
points(sp_locations)
points(sp_locations@coords[!i, 1],
sp_locations@coords[!i, 2],
col = "red", cex = 3, pch = 3
)
Now we can spot rejected points very easily as they are marked with a
red cross. If, however, we want to get rid of rejected points, we have
to apply a small trick: Applying the as.ppp() function on
the ppp object removes rejected points.
ppp_locations <- as.ppp(ppp_locations)
plot(ppp_locations)
No more warnings, great!
We continue by constructing ppp objects for the
berlin_random data set and the berlin_regular
data set. Since both data sets are data.frame objects we
simply index their columns with the \(\$\)-notation.
Exercise: Construct a
pppobject for theberlin_randomdata set and one forthe berlin_regulardata set, assign a unit for the coordinates and remove rejected points.
## Your code here...
ppp_random <- NULL
ppp_regular <- NULLShow code
## create ppp object
ppp_random <- ppp(
x = berlin_random$x,
y = berlin_random$y,
window = win
)## Warning: 184 points were rejected as lying outside the specified windowShow code
## assign unit
unitname(ppp_random) <- coordinate_unitsShow code
## removal of rejected points
ppp_random <- as.ppp(ppp_random)
plot(ppp_random)
Show code
## create ppp object
ppp_regular <- ppp(
x = berlin_regular$x,
y = berlin_regular$y,
window = win
)## Warning: 70 points were rejected as lying outside the specified windowShow code
## assign unit
unitname(ppp_regular) <- coordinate_unitsShow code
## removal of rejected points
ppp_regular <- as.ppp(ppp_regular)
plot(ppp_regular)
Finally, we save all three data sets for further usage.
save(
file = "ppp_data_berlin.RData",
ppp_locations, ppp_random, ppp_regular, berlin_locations
)Citation
The E-Learning project SOGA-R was developed at the Department of Earth Sciences by Kai Hartmann, Joachim Krois and Annette Rudolph. You can reach us via mail by soga[at]zedat.fu-berlin.de.
You may use this project freely under the Creative Commons Attribution-ShareAlike 4.0 International License.
Please cite as follow: Hartmann, K., Krois, J., Rudolph, A. (2023): Statistics and Geodata Analysis using R (SOGA-R). Department of Earth Sciences, Freie Universitaet Berlin.