Commit a3ea7bcb authored by Weigert, Andreas's avatar Weigert, Andreas
Browse files
parents a2d9bbf1 5a7b9673
.DS_Store
.Rproj.user .Rproj.user
.Rhistory .Rhistory
.RData .RData
......
---
title: "Practical tasks working with geographic and time-related data (EESYS-BIA-M, L05)"
output: html_notebook
editor_options:
chunk_output_type: inline
---
# Practicals for geographic data
```{r install packages}
install.packages("sp", "rgdal", "osmar", "lubridate")
```
Germany publishes geographic data as ESRI shapefiles on the Website of the [Bundesamt für Kartographie](http://www.geodatenzentrum.de/geodaten/gdz_rahmen.gdz_div?gdz_spr=deu&gdz_akt_zeile=5&gdz_anz_zeile=1&gdz_user_id=0)
For illusrration, the political borders (Verwaltungsgrenzen) are used on a scale of 1:250 000 (VG250)
There are two datasets available:
* VG250 Kompakt: contains the lines and general data (>120 MB)
* VG250 Ebenen: contains the detailed data (>70MB)
Selected layers are available in the `/data/` folder of this repository.
```{r load German political borders}
library(sp)
library(rgdal)
#read the shapefile
ger.shp <- readOGR("../data/DE_Admin_Geodata/VG250_Kompakt/", "VG250_L") #lines
#the attribute AGZ contains the type of the boundary (AGZ = 'Art der Grenze')
table(ger.shp$AGZ)
#Plot lines
plot(ger.shp[ger.shp$AGZ==1, ]) #German main border
lines(ger.shp[ger.shp$AGZ==2, ], col=2) #German state borders
lines(ger.shp[ger.shp$AGZ==3, ], col=3) #German Regierungsbezirk
lines(ger.shp[ger.shp$AGZ==4, ], col=4) #Kreise
```
```{r load German states as polygons}
#read teh shapefile
ger.shp <- readOGR("../data/DE_Admin_Geodata/VG250_Ebenen/", "VG250_LAN")
#basic plot of the data - looks similar to previous plot
plot(ger.shp)
axis(1)
axis(2)
#Here, we have polygons that we can colorate
plot(ger.shp, col=1:23)
#see the data stored in the shapefile
View(ger.shp@data)
```
As an example, the housing statistics from the German Federal Ministry of Statistics is used
Wohngebäude, Wohnungen, Wohnfläche: Bundesländer, Stichtag, Anzahl der Wohnungen
By this example, the `readr` package is illustrated
```{r visualize some data}
library(readr)
housing_statistic_DE <- read_csv2("../data/housing_statistic_germany_31231-0003.csv",
col_types = cols(
date = col_number(),
state = col_character(),
`1_app_buildings` = col_integer(),
`1_app_dwellings` = col_integer(),
`1_app_livingareas` = col_integer(),
`2_app_buildings` = col_integer(),
`2_app_dwellings` = col_integer(),
`2_app_livingareas` = col_integer(),
`3_more_app_buildings` = col_integer(),
`3_more_app_dwellings` = col_integer(),
`3_more_app_livingarea` = col_integer(),
dormitory_buildings = col_integer(),
dormitory_dwellings = col_integer(),
dormitory_livingarea = col_integer(),
total_buildings = col_integer(),
total_dwellings = col_integer(),
total_livingarea = col_integer()
),
col_names = c("date", "state",
"1_app_buildings", "1_app_dwellings", "1_app_livingareas",
"2_app_buildings", "2_app_dwellings", "2_app_livingareas",
"3_more_app_buildings", "3_more_app_dwellings", "3_more_app_livingarea",
"dormitory_buildings", "dormitory_dwellings", "dormitory_livingarea",
"total_buildings", "total_dwellings", "total_livingarea"),
skip = 8, n_max = 16, locale = locale(date_names = "de",
encoding = "latin1",
decimal_mark = ",",
grouping_mark = "."))
View(housing_statistic_DE)
#calculate the relative distribution of total living area
housing_statistic_DE$total_livingarea_rel <- housing_statistic_DE$total_livingarea / max(housing_statistic_DE$total_livingarea)
#the HSV function definex colors by means of hue, saturation and darkness value
color_grades <- hsv(h=0.6, s=housing_statistic_DE$total_livingarea_rel, v = 1)
pos_states_map_statistics <- match(ger.shp$GEN, housing_statistic_DE$state)
plot(ger.shp, col=color_grades[pos_states_map_statistics],
main="Total living space area in the German states")
axis(1)
axis(2)
```
```{r plot coordinaes on the map}
# the WGS84 coordinates of Bamberg are
bamberg <- SpatialPoints(list(x=10.891667, y=49.891667))
#display the projection and coordinate system of the spatial data
proj4string(ger.shp)
proj4string(bamberg)
# transform the german UTM coordinates to WGS84 coordinates
ger.shp.wgs84 <- spTransform(ger.shp, CRS("+proj=longlat +datum=WGS84"))
plot(ger.shp.wgs84)
axis(1)
axis(2)
points(bamberg, col=2)
```
Here is an example on how to access OpenStreetMap
**Please** do not query OSM with more than 1 request per second!!!
The API is not made for large read queries
There is a special [Overpass API](https://wiki.openstreetmap.org/wiki/Overpass_API)
to retrieve much data from the database.
```{r accessing OpenStreetMap data}
library(osmar)
erba_data <- get_osm(x = center_bbox(center_lon = 10.86832,
center_lat = 49.90300,
width=300,
height=300),
source = osmsource_api(url = "https://api.openstreetmap.org/api/0.6/"))
#plot the data
plot(as_sp(erba_data, what = "lines"))
points(as_sp(erba_data, what = "points"), pch=16, col=2, cex=0.5)
axis(1)
axis(2)
#display the data stored for each object
View(erba_data$nodes$attrs)
#display the tags stored for each object (contain the geographic meaning)
View(erba_data$nodes$attrs)
```
# Practicals for the date and time handling
```{r parsing data}
library(lubridate)
#for the case you get a warning with "unknown timezone"
Sys.setenv(TZ="CET")
dmy(13112018)
ymd("20181113")
mdy(c("November 13th, 2018", "11/ 13/ 2018"))
dmy_hm("13.11.2018 11:15")
dmy_hms("13.11.2018 11:15:01", tz = "CET")
today()
now()
ymd(c("2010-10-10", "bananas"))
```
```{r getting and setting components}
t <- dmy_hms("13.11.2018 11:15:01", tz = "CET")
second(t)
minute(t)
day(t)
week(t)
month(t)
month(t, label = T)
wday(t)
wday(t, label = T)
year(t) <- 2017
print(t)
#override time zone
tz(t) <- "EST"
t
#converts the time zone
with_tz(t, "CET")
```
```{r date time aritmetics}
t <- c("13.11.2018 11:15:01",
"12.11.2018 23:16:09",
"09.11.2018 07:01:01",
"08.11.2018 10:19:10")
t <- dmy_hms(t)
mean(t)
floor_date(t, unit = "week")
```
There are **three types of time spans**:
* durations, which represent an exact number of seconds.
* periods, which represent human units like weeks and months.
* intervals, which represent a starting and ending point.
`difftime` objects are hard to work with as they store the difference in seconds, minutes, ...
```{r time spans}
#how old is Hannah?
h_age <- today() - dmy(21031985)
h_age
class(h_age) #duration stored as a difftim-object
h_age_dur <- as.duration(h_age) #lubridate offers a duration object
h_age_dur
#durations can also be created
dweeks(3)
dyears(2)
#mind the fact:
dmy_hms("27.10.2018 10:00:00", tz="CET")
dmy_hms("27.10.2018 10:00:00", tz="CET")+ddays(1) #time shift for DST
#to solve this, lubridate offers periods
dmy_hms("27.10.2018 10:00:00", tz="CET")+days(1)
# an example with leap years
ymd("2016-01-01") + dyears(1)
ymd("2016-01-01") + years(1)
```
```{r intervals}
next_year <- today() + years(1)
#lubridate brings a new operation for time intervals
today() %--% next_year
(today() %--% next_year) / ddays(1)
#there is also an operation to check how many periods fall into an interval
(today() %--% next_year) %/% dweeks(1)
```
PROJCS["ETRS_1989_UTM_Zone_32N",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
PROJCS["ETRS_1989_UTM_Zone_32N",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
UTF-8
\ No newline at end of file
PROJCS["ETRS_1989_UTM_Zone_32N",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment