Analyzing environmental sensor data from openSenseMap.org in R
This package provides data ingestion functions for almost any data stored on the open data platform for environemental sensordata https://opensensemap.org. Its main goals are to provide means for:
- big data analysis of the measurements stored on the platform
- sensor metadata analysis (sensor counts, spatial distribution, temporal trends)
Exploring the dataset
Before we look at actual observations, lets get a grasp of the openSenseMap datasets’ structure.
library(magrittr)
library(opensensmapr)
all_sensors = osem_boxes()summary(all_sensors)## boxes total: 704
##
## boxes by exposure:
## indoor outdoor unknown
## 127 556 21
##
## boxes by model:
## custom homeEthernet homeEthernetFeinstaub
## 209 78 8
## homeWifi homeWifiFeinstaub luftdaten_sds011
## 106 34 22
## luftdaten_sds011_bme280 luftdaten_sds011_bmp180 luftdaten_sds011_dht11
## 41 3 14
## luftdaten_sds011_dht22
## 189
##
## $last_measurement_within
## 1h 1d 30d 365d never
## 0 0 0 0 704
##
## oldest box: 2014-05-28 15:36:14 (CALIMERO)
## newest box: 2017-08-24 17:38:44 (Burgweinting)
##
## sensors per box:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 4.000 5.000 4.605 5.000 17.000This gives a good overview already: As of writing this, there are more than 700 sensor stations, of which ~50% are currently running. Most of them are placed outdoors and have around 5 sensors each. The oldest station is from May 2014, while the latest station was registered a couple of minutes ago.
Another feature of interest is the spatial distribution of the boxes: plot() can help us out here. This function requires a bunch of optional dependencies though.
if (!require('maps')) install.packages('maps')
if (!require('maptools')) install.packages('maptools')
if (!require('rgeos')) install.packages('rgeos')
plot(all_sensors)
It seems we have to reduce our area of interest to Germany.
But what do these sensor stations actually measure? Lets find out. osem_phenomena() gives us a named list of of the counts of each observed phenomenon for the given set of sensor stations:
phenoms = osem_phenomena(all_sensors)
str(phenoms)## List of 189
## $ Temperatur : int 647
## $ rel. Luftfeuchte : int 534
## $ Luftdruck : int 368
## $ PM10 : int 347
## $ PM2.5 : int 347
## $ UV-Intensität : int 255
## $ Beleuchtungsstärke : int 251
## $ Luftfeuchtigkeit : int 83
## $ Schall : int 26
## $ Helligkeit : int 20
## $ Licht : int 20
## $ UV : int 15
## $ Humidity : int 12
## $ Temperature : int 11
## $ Anderer : int 10
## $ Ilmanpaine : int 9
## $ Lämpötila : int 9
## $ Licht (digital) : int 9
## $ Valonmäärä : int 8
## $ Windgeschwindigkeit : int 8
## $ Kosteus : int 7
## $ Luftfeuchte : int 7
## $ Lautstärke : int 6
## $ Signal : int 6
## $ UV-säteily : int 6
## $ Wind speed : int 5
## $ Pressure : int 4
## $ temperature : int 4
## $ Windrichtung : int 4
## $ DS18B20_Probe01 : int 3
## $ DS18B20_Probe02 : int 3
## $ DS18B20_Probe03 : int 3
## $ DS18B20_Probe04 : int 3
## $ DS18B20_Probe05 : int 3
## $ Light : int 3
## $ Niederschlag : int 3
## $ UV Index : int 3
## $ UV-Säteily : int 3
## $ UV-Strahlung : int 3
## $ C2H5OH : int 2
## $ CO : int 2
## $ CPU-Temp : int 2
## $ Feinstaub : int 2
## $ Feinstaub PM10 : int 2
## $ Feinstaub PM2,5 : int 2
## $ H2 : int 2
## $ humidity : int 2
## $ Ilmankosteus : int 2
## $ NH3 : int 2
## $ NO2 : int 2
## $ Regen : int 2
## $ rel. Luftfeuchtigkeit : int 2
## $ Relative Humidity : int 2
## $ Temperatur BMP280 : int 2
## $ Temperatur DHT22 : int 2
## $ Temperatur HDC1008 : int 2
## $ TemperaturBME : int 2
## $ test : int 2
## $ UV-Index : int 2
## $ Wassertemperatur : int 2
## $ Wifi-Stärke : int 2
## $ Windböen : int 2
## $ Wolkenbedeckung : int 2
## $ Air Preassure : int 1
## $ Air pressure : int 1
## $ Air Temperature : int 1
## $ Akkuspannung Terrasse : int 1
## $ Akkuspannung Unten Eingang : int 1
## $ Attendance : int 1
## $ Batterie : int 1
## $ Batterieladung : int 1
## $ Battery : int 1
## $ Beleuchtungsstaerke : int 1
## $ Beleuchtungsstärke des sichtbaren Lichts: int 1
## $ Bodenfeuchte : int 1
## $ Bodentemperatur : int 1
## $ C3H8 : int 1
## $ C4H10 : int 1
## $ CH4 : int 1
## $ CO2 : int 1
## $ CO2-Konzentration : int 1
## $ Dämmerung : int 1
## $ dT : int 1
## $ Dust Sensor : int 1
## $ Dust_Concentration : int 1
## $ Eingangsspannung : int 1
## $ Feinstaub P10 : int 1
## $ Feinstaub P2.5 : int 1
## $ Feinstaubgehalt PM10 : int 1
## $ Feinstaubgehalt PM2.5 : int 1
## $ Feinstaubkonzentration : int 1
## $ Feuchte : int 1
## $ Feuchtigkeit : int 1
## $ filedData : int 1
## $ H2, LPG, CH4, CO, Alcohol : int 1
## $ Höhe : int 1
## $ Illuminance : int 1
## $ Infrared light : int 1
## $ Intensität der ultravioletten Strahlung : int 1
## [list output truncated]Thats quite some noise there, with many phenomena being measured by a single sensor only, or many duplicated phenomena due to slightly different spellings. We should clean that up, but for now let’s just filter out the noise and find those phenomena with high sensor numbers:
phenoms[phenoms > 20]## $Temperatur
## [1] 647
##
## $`rel. Luftfeuchte`
## [1] 534
##
## $Luftdruck
## [1] 368
##
## $PM10
## [1] 347
##
## $PM2.5
## [1] 347
##
## $`UV-Intensität`
## [1] 255
##
## $Beleuchtungsstärke
## [1] 251
##
## $Luftfeuchtigkeit
## [1] 83
##
## $Schall
## [1] 26Alright, temperature it is! Fine particulate matter (PM2.5) seems to be more interesting to analyze though. We should check how many sensor stations provide useful data: We want only those boxes with a PM2.5 sensor, that are placed outdoors and are currently submitting measurements:
pm25_sensors = osem_boxes(
exposure = 'outdoor',
date = Sys.time(), # ±4 hours
phenomenon = 'PM2.5'
)summary(pm25_sensors)## boxes total: 240
##
## boxes by exposure:
## outdoor
## 240
##
## boxes by model:
## custom homeEthernetFeinstaub homeWifi
## 18 4 5
## homeWifiFeinstaub luftdaten_sds011 luftdaten_sds011_bme280
## 12 14 29
## luftdaten_sds011_bmp180 luftdaten_sds011_dht11 luftdaten_sds011_dht22
## 1 11 146
##
## $last_measurement_within
## 1h 1d 30d 365d never
## 0 0 0 0 240
##
## oldest box: 2016-09-11 08:17:17 (Balkon Gasselstiege)
## newest box: 2017-08-24 17:38:44 (Burgweinting)
##
## sensors per box:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.000 4.000 4.000 4.275 4.000 10.000plot(pm25_sensors)
Thats still more than 200 measuring stations, we can work with that.
Analyzing sensor data
Having analyzed the available data sources, let’s finally get some measurements. We could call osem_measurements(pm25_sensors) now, however we are focussing on a restricted area of interest, the city of Berlin. Luckily we can get the measurements filtered by a bounding box:
library(sf)
library(units)
library(lubridate)
library(dplyr)
# construct a bounding box: 12 kilometers around Berlin
berlin = st_point(c(13.4034, 52.5120)) %>%
st_sfc(crs = 4326) %>%
st_transform(3857) %>% # allow setting a buffer in meters
st_buffer(set_units(12, km)) %>%
st_transform(4326) %>% # the opensensemap expects WGS 84
st_bbox()pm25 = osem_measurements(
berlin,
phenomenon = 'PM2.5',
from = now() - days(20), # defaults to 2 days
to = now()
)
plot(pm25)
Now we can get started with actual spatiotemporal data analysis. First, lets mask the seemingly uncalibrated sensors:
outliers = filter(pm25, value > 100)$sensorId
bad_sensors = outliers[, drop = T] %>% levels()
pm25 = mutate(pm25, invalid = sensorId %in% bad_sensors)Then plot the measuring locations, flagging the outliers:
st_as_sf(pm25) %>% st_geometry() %>% plot(col = factor(pm25$invalid), axes = T)
Removing these sensors yields a nicer time series plot:
pm25 %>% filter(invalid == FALSE) %>% plot()
Further analysis: comparison with LANUV data TODO