Introduction - Class: SpatialPoint 2

This exercise helps to get a good feeling on how to handle
- shape files
- relief files
- transforming coordinates into another reference system

Three different data are being processed here:
- Municipality geometries: The geometries do not show the political borders of Swiss municipalities, but the so-called “productive” area, i.e., larger lakes and other “unproductive” areas such as mountains are excluded.
- Earth quakes data from the year 1984, filtered out the data for Switzerland.
- and a freely available GeoTIFF (relief)

Important: I am not answering any relevant questions here. This exercise shall merely help me to get a good grasp on how to handle these different file types.

remove(list = ls(all.names = TRUE))
cat("\014")


require(dplyr)
require(broom)
require(rgdal)
require(raster)
require(maptools)
# require(xtable) # for a few nice tables
require(ggplot2)
require(gridExtra)
# data <- read.csv("input/avg_age_15.csv", stringsAsFactors = F) 
# dim(data)
# head(data)

Read in geodata

Details about CRS are given in my previous exercise about SpatialPoints; also about who to handle the this class SpatialPoint.
- readOGR reads in the shapefile (.shp).
- Also important is to note, that the function tidy needs to be used, instaead of fortify*. For this, I needed two packages: broom and a dependency of it maptools. tidy makes the shapefile ggplot2 compatible.

gde_15 <- readOGR("input/geodata/gde-1-1-15.shp", layer = "gde-1-1-15")
#> OGR data source with driver: ESRI Shapefile 
#> Source: "input/geodata/gde-1-1-15.shp", layer: "gde-1-1-15"
#> with 2324 features
#> It has 2 fields
#> Integer64 fields read as strings:  BFS_ID

# set crs to ch1903/lv03, just to make sure  (EPSG:21781)
crs(gde_15) <- "+proj=somerc +lat_0=46.95240555555556 
+lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 
+ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"

# fortify, i.e., make ggplot2-compatible
#' new way to 'fortify' with broom

# library(broom)
map_data_fortified <- tidy(gde_15, region = "BFS_ID")  %>% 
          mutate(id = as.numeric(id))
#> Warning: package 'bindrcpp' was built under R version 3.3.3

Let’s have a glimpse on how the data looks like.
Keep i nmind, the “fortified” data is a data frame, not a SpatialPoint class (not yet).

dim(map_data_fortified)
#> [1] 29040     7
class(map_data_fortified)
#> [1] "data.frame"
head(map_data_fortified)
#>       long    lat order  hole piece group id
#> 1 678565.1 238072     1 FALSE     1   1.1  1
#> 2 679757.1 238500     2 FALSE     1   1.1  1
#> 3 680802.1 237767     3 FALSE     1   1.1  1
#> 4 680610.1 236347     4 FALSE     1   1.1  1
#> 5 680544.1 235029     5 FALSE     1   1.1  1
#> 6 680103.1 235061     6 FALSE     1   1.1  1

For now, to keep the exercise simple, let’s give each different canton a color. In the data frame, each municipality has an id and each canton has a defined range of id’s. This you can look up in Wikipedia, for instance.


map_data <- map_data_fortified # %>% left_join(data, by = c("id" = "bfs_id"))
dim(map_data)
#> [1] 29040     7
head(map_data)
#>       long    lat order  hole piece group id
#> 1 678565.1 238072     1 FALSE     1   1.1  1
#> 2 679757.1 238500     2 FALSE     1   1.1  1
#> 3 680802.1 237767     3 FALSE     1   1.1  1
#> 4 680610.1 236347     4 FALSE     1   1.1  1
#> 5 680544.1 235029     5 FALSE     1   1.1  1
#> 6 680103.1 235061     6 FALSE     1   1.1  1

map_data <- within(map_data, {
        canton <- NA
        canton[id >= 0001 & id <= 0298]         <- "Zurich" # original: 0261
        canton[id >= 0301 & id <= 0996]         <- "Bern"
        canton[id >= 1001 & id <= 1151]         <- "Luzern" # Original: 1150
        canton[id >= 1201 & id <= 1220]         <- "Uri"
        canton[id >= 1301 & id <= 1375]         <- "Schwyz"
        canton[id >= 1401 & id <= 1407]         <- "Obwalden"
        canton[id >= 1501 & id <= 1511]         <- "Nidwalden"
        canton[id >= 1601 & id <= 1632]         <- "Glarus" # original: 1629
        canton[id >= 1701 & id <= 1711]         <- "Zug"
        canton[id >= 2001 & id <= 2338]         <- "Freiburg" # original: 2336
        canton[id >= 2401 & id <= 2622]         <- "Solothurn"
        canton[id >= 2701 & id <= 2703]         <- "Basel-Stadt"
        canton[id >= 2761 & id <= 2895]         <- "Basel-Landschaft"
        canton[id >= 2901 & id <= 2974]         <- "Schaffhausen"
        canton[id >= 3001 & id <= 3038]         <- "Appenzell Ausserrhoden"
        canton[id >= 3101 & id <= 3111]         <- "Appenzell Innerrhoden"
        canton[id >= 3201 & id <= 3444]         <- "St. Gallen"
        canton[id >= 3501 & id <= 3987]         <- "Graubuenden"
        canton[id >= 4001 & id <= 4323]         <- "Aargau"
        canton[id >= 4401 & id <= 4951]         <- "Thurgau"
        canton[id >= 5001 & id <= 5398]         <- "Tessin" # original: 5322
        canton[id >= 5401 & id <= 5939]         <- "Waadt" # original: 5939 
        canton[id >= 6001 & id <= 6300]         <- "Wallis"
        canton[id >= 6401 & id <= 6512]         <- "Neuenburg" # original 6511
        canton[id >= 6601 & id <= 6645]         <- "Genf"
        canton[id >= 6701 & id <= 6810]         <- "Jura" # original 6806
        
})



# Interestingly, there were a few municipalities, which were out of the defined range. 
# Why?
# I don't know yet.
#
# on my to-do-list 

sum(is.na(map_data$canton))
#> [1] 0
# map_data_NA <- map_data[is.na(map_data$canton), ]
# table(map_data_NA$id)

There were a couple of out-of-the-defined-range municipalities. I don’t yet know why. However, you can see in the code above, that I slightly modified the id-range, to incorporate all ids.

First Plot

Very straight forward: how does it look like.

p1 <- ggplot() +
        geom_polygon(data = map_data, aes(fill = canton, 
                                          x = long, 
                                          y = lat, 
                                          group = group)) + # 
        geom_path(data = map_data, aes(long, lat, group = group))


p1

Importing the relief

WORK IN PROGRESS
- read in the swiss relief
- relief has geo-data
- put this relief beneath the color plot

Earth Quake Exercise

You remember the exercise 1 of mine during which we took the very same data file of earthquakes from the year 1984. Here, we repeate the very same procedure.

data <- read.table("./010_data/1984.txt", sep = "|", header = TRUE, comment.char = "", quote = "")
data <- filter(data, EventLocationName == "SWITZERLAND")
dim(data)
#> [1] 61 13
head(data)
#>   X.EventID                Time Latitude Longitude Depth.km Author Catalog
#> 1   2971993 1984-01-11T14:11:57  47.3872    8.7689     10.0    ISC     ISC
#> 2   2977978 1984-01-19T04:01:54  47.5792    7.6617      6.0    ISC     ISC
#> 3   2978110 1984-01-19T08:36:45  46.1723    9.2092     13.0    ISC     ISC
#> 4   2978919 1984-01-20T02:57:11  46.5968    8.4459      6.1    ISC     ISC
#> 5   2983336 1984-01-26T00:01:36  47.3000    8.0000       NA    LDG     ISC
#> 6   2984861 1984-01-28T02:35:35  47.5136    8.2287      0.0    ISC     ISC
#>   Contributor ContributorID MagType Magnitude MagAuthor EventLocationName
#> 1         ISC        560836      mL       3.7       LDG       SWITZERLAND
#> 2         ISC        561325      mL       1.9       ZUR       SWITZERLAND
#> 3         ISC        561338      mL       2.2       ZUR       SWITZERLAND
#> 4         ISC        561401      mL       2.1       ZUR       SWITZERLAND
#> 5         ISC        561760      mL       2.2       LDG       SWITZERLAND
#> 6         ISC        561897      mL       2.4       LDG       SWITZERLAND

Defining the new CRS:

CRS.new <- CRS("+proj=somerc +lat_0=46.9524056 +lon_0=7.43958333 +ellps=bessel +x_0=600000 +y_0=200000 +towgs84=674.374,15.056,405.346 +units=m +k_0=1 +no_defs")

side note: the following is another CRS definition I found on the internet. It clearly has another x_0 and y_0. I wonder, which CRS this is.
CRS.new <- CRS(“+proj=somerc +lat_0=46.9524056 +lon_0=7.43958333 +ellps=bessel +x_0=2600000 +y_0=1200000 +towgs84=674.374,15.056,405.346 +units=m +k_0=1 +no_defs”)

source: [http://gis.stackexchange.com/questions/45263/converting-geographic-coordinate-system-in-r]

Here the very same routine as in the previous exercise:

coordinates(data) <-  c("Longitude", "Latitude")  # change into SpatialPoint
proj4string(data) <- CRS("+init=epsg:4326") # point of reference = WGS 84 (World Geodetic System 1984)
head(data) # Let's see how it looks like
#>   X.EventID                Time Depth.km Author Catalog Contributor
#> 1   2971993 1984-01-11T14:11:57     10.0    ISC     ISC         ISC
#> 2   2977978 1984-01-19T04:01:54      6.0    ISC     ISC         ISC
#> 3   2978110 1984-01-19T08:36:45     13.0    ISC     ISC         ISC
#> 4   2978919 1984-01-20T02:57:11      6.1    ISC     ISC         ISC
#> 5   2983336 1984-01-26T00:01:36       NA    LDG     ISC         ISC
#> 6   2984861 1984-01-28T02:35:35      0.0    ISC     ISC         ISC
#>   ContributorID MagType Magnitude MagAuthor EventLocationName
#> 1        560836      mL       3.7       LDG       SWITZERLAND
#> 2        561325      mL       1.9       ZUR       SWITZERLAND
#> 3        561338      mL       2.2       ZUR       SWITZERLAND
#> 4        561401      mL       2.1       ZUR       SWITZERLAND
#> 5        561760      mL       2.2       LDG       SWITZERLAND
#> 6        561897      mL       2.4       LDG       SWITZERLAND
data.ch1903 <- spTransform(data, CRS.new)
head(data.ch1903)
#>   X.EventID                Time Depth.km Author Catalog Contributor
#> 1   2971993 1984-01-11T14:11:57     10.0    ISC     ISC         ISC
#> 2   2977978 1984-01-19T04:01:54      6.0    ISC     ISC         ISC
#> 3   2978110 1984-01-19T08:36:45     13.0    ISC     ISC         ISC
#> 4   2978919 1984-01-20T02:57:11      6.1    ISC     ISC         ISC
#> 5   2983336 1984-01-26T00:01:36       NA    LDG     ISC         ISC
#> 6   2984861 1984-01-28T02:35:35      0.0    ISC     ISC         ISC
#>   ContributorID MagType Magnitude MagAuthor EventLocationName
#> 1        560836      mL       3.7       LDG       SWITZERLAND
#> 2        561325      mL       1.9       ZUR       SWITZERLAND
#> 3        561338      mL       2.2       ZUR       SWITZERLAND
#> 4        561401      mL       2.1       ZUR       SWITZERLAND
#> 5        561760      mL       2.2       LDG       SWITZERLAND
#> 6        561897      mL       2.4       LDG       SWITZERLAND

Let’s see the data in the two reference system:

par(mfrow = c(1, 2))
# plot.default(x, y, main = "Raw data - x and y", cex.axis = .95)
plot(data$Longitude, data$Latitude, axes = TRUE, main = "Original lat-lon - d", cex.axis = .95)
plot(data.ch1903$Longitude, data.ch1903$Latitude, axes = TRUE, main = "Projected - d.ch1903", cex.axis = .95)

# unclass(d.ch1903)

Very rough now: let’s plot the earth quake locations on the swiss map - just for fun:

p2 <- p1 + geom_point(aes(data.ch1903$Longitude, data.ch1903$Latitude, 
                          size = data.ch1903$Magnitude)) +
        theme(legend.position="none")
p2

One can ask, if I really made the CRS-transformation correctly. Fair question. For that, let’s go back and have a look at the original longitudinal and latitudinal data in the data frame data.
Below, you see the plot of the world, zooming in with coord_map and plotting the earth quake data from data. Bear in mind, data has the unaltered coordinates in it. Only the data frame data.ch1903 has the projected ones in it.

ggP <- ggplot() + 
        borders("world", colour="black", fill="#FFF8DC") + 
        xlab("") + #remove the x- and y-axis labels# create a layer of borders
        ylab("") +
        coord_map(xlim = c(5.5, 11), ylim = c(45.5, 48)) + #zoom EU  
        geom_point(aes(data$Longitude, data$Latitude)
                   , shape = 21
                   # , data = data
                   , color = "red" #reddish color
                   , fill="salmon1"  #good ones: yellow1 , sienna1 , palevioletred2 , salmon1
                   , alpha = 0.3
        )
#> Warning: package 'maps' was built under R version 3.3.3

print(ggP)

Two put both plots next to each other, I needed the package gridExtra, because ggplot2 isn’t well equipped in this matter.
One sees, that the data points are the same ones, independent on the coordinate system. So we conclude, that we transformed the data correctly

gridExtra::grid.arrange(p2, ggP, nrow = 1)

Importing the relief

WORK IN PROGRESS
We saw that the transformation was done correctly. We can now get rid of the second plot (yellowish), which only served to have a comparison to the original data values.
First, let’s load the GeoTIF and convert it to a spatial pixel data frame (spdf); and then extract its content with using as.data.frame:

relief <- raster("input/geodata/02-relief-georef-clipped-resampled.tif")
relief_spdf <- as(relief, "SpatialPixelsDataFrame")
relief <- as.data.frame(relief_spdf) %>% 
        rename(value = `X02.relief.georef.clipped.resampled`)
dim(relief)
#> [1] 1135311       3
head(relief)
#>   value        x        y
#> 1   230 687087.7 294707.0
#> 2   230 687276.4 294707.0
#> 3   230 687465.2 294707.0
#> 4   214 685577.7 294518.2
#> 5   214 685766.4 294518.2
#> 6   215 685955.2 294518.2

To incorporate the relief into the Swiss map, we take up the map p1 - removing the legend and adding the relief via the geom_polygon command:

p1 <- ggplot() +
        # relief is the background
        geom_raster(data = relief, aes(x = x, 
                                       y = y, 
                                       alpha = value)) +
        # plotting the cantons:
        geom_polygon(data = map_data, aes(fill = canton, 
                                          x = long, 
                                          y = lat, 
                                          group = group)) + # 
        geom_path(data = map_data, aes(long, lat, group = group)) +
        theme(legend.position="none")


p1

It seems quite dark. Let’s take the same code and add some alpha value correction for the relief:

p1 <- ggplot() +
        # relief is the background
        geom_raster(data = relief, aes(x = x, 
                                       y = y, 
                                       alpha = value)) +
        # alpha value correction
        scale_alpha(range = c(0.8, 0), guide = F)  + 
        # plotting the cantons:
        geom_polygon(data = map_data, aes(fill = canton, 
                                          x = long, 
                                          y = lat, 
                                          group = group)) + # 
        geom_path(data = map_data, aes(long, lat, group = group)) +
        theme(legend.position="none")


p1

Conclusion

Of course, we can take this plot even further, much further.
However, for the purpose of this exercise, we stop here.
- We transformed the coordinates successfully,
- The Swiss shape files can be seen in the plot, according to the municipalities.
- And we added a GeoTIF-relief underneath.
It looks neat.