Introduction - Class: SpatialPoint

This exercise shall give me some insights on how to get from “normal” longitudes and latitudes to another coordinate reference system (in this case, the swiss grid).

  • The package rgdal gives us access to projection/transformation operations from the PROJ.4 library. Both GDAL raster and OGR vector map data can be imported into R, and GDAL raster data and OGR vector data exported.
  • And as usual, dplyr for data manipulation
require(rgdal)
require(dplyr)

Often, I have longitudes and latitudes given in a similar data frame as the one below:

x <- c(7.173500, 7.172540, 7.171636, 7.180180, 7.178070, 7.177229, 7.175240, 7.181409, 7.179299)
y <- c(45.86880, 45.86887, 45.86924, 45.87158, 45.87014, 45.86923, 45.86808, 45.87177, 45.87020)

(d <- data.frame(lon=x, lat=y))
#>        lon      lat
#> 1 7.173500 45.86880
#> 2 7.172540 45.86887
#> 3 7.171636 45.86924
#> 4 7.180180 45.87158
#> 5 7.178070 45.87014
#> 6 7.177229 45.86923
#> 7 7.175240 45.86808
#> 8 7.181409 45.87177
#> 9 7.179299 45.87020

To changing the data frame d into a SpatialPoints-class, merely define coordinates to it, and voilà. Additioanlly, the coordinate reference system for these cooridantes is essential to have.

coordinates(d) <- c("lon", "lat")  # change into SpatialPoint
proj4string(d) <- CRS("+init=epsg:4326") # point of reference = WGS 84 (World Geodetic System 1984)
d # Let's see how it looks like
#> SpatialPoints:
#>            lon      lat
#>  [1,] 7.173500 45.86880
#>  [2,] 7.172540 45.86887
#>  [3,] 7.171636 45.86924
#>  [4,] 7.180180 45.87158
#>  [5,] 7.178070 45.87014
#>  [6,] 7.177229 45.86923
#>  [7,] 7.175240 45.86808
#>  [8,] 7.181409 45.87177
#>  [9,] 7.179299 45.87020
#> Coordinate Reference System (CRS) arguments: +init=epsg:4326
#> +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

4326 is just the EPSG identifier of WGS84 - which is the world geodetic system defined in 1984. WGS84 comprises a standard coordinate frame for the Earth, a datum/reference ellipsoid for raw altitude data, and a gravitational equipotential surface (the geoid) that defines the nominal sea level.

The new CRS shall be the Swiss Grid, which is defined by
- fundamental point the old observatory of Bern (46°57′3.9″N 7°26′19.1″E (WGS84)) - the coordinates of this point are 600’000 m E / 200’000 m N
- The 0 / 0 coordinate is located near Bordeaux, France
- Though E coordinate is denoted as y and N coordinate x, E coordinate is the first axis of this Cartesian system, namely a point is denoted as (y, x). As an example, the observatory is (600’000, 200’000)
- All coordinates are always positive, since Switzerland is located in the 1st quadrant of the coordinate system.
- Furthermore, the whole area of Switzerland is located below the y=x line of the coordinate system. Thus, all E-coordinates are always bigger than N-coordinates.
- The map projection used is Oblique Mercator on an 1841 Bessel ellipsoid.

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 was the origianl CRS, which I found researching for it:
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]

Now let’s take the SpatialPoints ( = d) and transform into the new CRS:

(d.ch1903 <- spTransform(d, CRS.new))
#> SpatialPoints:
#>            lon      lat
#>  [1,] 579408.2 79721.15
#>  [2,] 579333.7 79729.18
#>  [3,] 579263.7 79770.55
#>  [4,] 579928.0 80028.46
#>  [5,] 579763.6 79868.92
#>  [6,] 579698.0 79767.97
#>  [7,] 579543.1 79640.65
#>  [8,] 580023.5 80049.26
#>  [9,] 579859.1 79875.27
#> Coordinate Reference System (CRS) arguments: +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

Plots

What happened graphically? Let’s compare the original data points with the new ones:
- original data points: vecotr x and y
- we transformed the data frame d into the WGS84 (carefull with the y-axis while comparing)
- and the final data frame d.ch1903

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

unclass(d.ch1903)
#> <S4 Type Object>
#> attr(,"bbox")
#>           min       max
#> lon 579263.65 580023.55
#> lat  79640.65  80049.26
#> attr(,"proj4string")
#> CRS arguments:
#>  +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 
#> attr(,"coords")
#>            lon      lat
#>  [1,] 579408.2 79721.15
#>  [2,] 579333.7 79729.18
#>  [3,] 579263.7 79770.55
#>  [4,] 579928.0 80028.46
#>  [5,] 579763.6 79868.92
#>  [6,] 579698.0 79767.97
#>  [7,] 579543.1 79640.65
#>  [8,] 580023.5 80049.26
#>  [9,] 579859.1 79875.27

Testing - with a dataset of earthquakes

In another project of mine, I was working on a earth quake data set. I took from there the data from 1984, given as a txt file and filtering out the quakes within Switzerland.
Read the data in - and have a look at it:

data <- read.table("./010_data/1984.txt", sep = "|", header = TRUE, comment.char = "", quote = "")
data <- filter(data, EventLocationName == "SWITZERLAND")
#> Warning: package 'bindrcpp' was built under R version 3.3.3
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

As before, define where the coordinates are in the data frame and define the CRS of them:

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
#>         coordinates X.EventID                Time Depth.km Author Catalog
#> 1 (8.7689, 47.3872)   2971993 1984-01-11T14:11:57     10.0    ISC     ISC
#> 2 (7.6617, 47.5792)   2977978 1984-01-19T04:01:54      6.0    ISC     ISC
#> 3 (9.2092, 46.1723)   2978110 1984-01-19T08:36:45     13.0    ISC     ISC
#> 4 (8.4459, 46.5968)   2978919 1984-01-20T02:57:11      6.1    ISC     ISC
#> 5         (8, 47.3)   2983336 1984-01-26T00:01:36       NA    LDG     ISC
#> 6 (8.2287, 47.5136)   2984861 1984-01-28T02:35:35      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

Same transformation - with new CRS:

data.ch1903 <- spTransform(data, CRS.new)
head(data.ch1903)
#>            coordinates X.EventID                Time Depth.km Author
#> 1 (700441.1, 249336.7)   2971993 1984-01-11T14:11:57     10.0    ISC
#> 2 (616783.4, 269857.2)   2977978 1984-01-19T04:01:54      6.0    ISC
#> 3 (736728.5, 114967.8)   2978110 1984-01-19T08:36:45     13.0    ISC
#> 4   (677182.6, 161110)   2978919 1984-01-20T02:57:11      6.1    ISC
#> 5   (642457, 238942.3)   2983336 1984-01-26T00:01:36       NA    LDG
#> 6 (659514.2, 262838.7)   2984861 1984-01-28T02:35:35      0.0    ISC
#>   Catalog Contributor ContributorID MagType Magnitude MagAuthor
#> 1     ISC         ISC        560836      mL       3.7       LDG
#> 2     ISC         ISC        561325      mL       1.9       ZUR
#> 3     ISC         ISC        561338      mL       2.2       ZUR
#> 4     ISC         ISC        561401      mL       2.1       ZUR
#> 5     ISC         ISC        561760      mL       2.2       LDG
#> 6     ISC         ISC        561897      mL       2.4       LDG
#>   EventLocationName
#> 1       SWITZERLAND
#> 2       SWITZERLAND
#> 3       SWITZERLAND
#> 4       SWITZERLAND
#> 5       SWITZERLAND
#> 6       SWITZERLAND

And a small plot, comparing before and after transforming

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)

Final Words

We will follow up this exercise in the next one: SpatialPoints 2