plot_data.R

gistfile1.R

library(dplyr)
library(tidyr)
library(maps)
library(rworldmap)
library(mapdata)
library(rgdal)
library(shape)

medians <- read.csv("wintershifts_median.csv", header = TRUE, sep = ";",
                    stringsAsFactors = FALSE)

# Määritellään LAEA koordinaattijärjestelmä
crs.laea <- CRS("+init=epsg:3035")

# Muuta spatiaaliseksi dataksi, CRS on latlon
sp_data <- SpatialPointsDataFrame(coords = select(medians, Long, Lat),
                                  data = select(medians, species, ID, year), 
                                  proj4string = CRS("+init=epsg:4326"))
# Projisoi LAEA koordinaatistoon
sp_data_laea <- spTransform(sp_data, crs.laea)
# Irrotetaan atribuuttidata uusilla koordinaateilla
sp_data_laea <- as.data.frame(sp_data_laea)

# Tarvitaan jottei as.double hukkaa tarkkuutta
options(digits = 16)
# Puljaa dataa
sp_data_laea_change <- sp_data_laea %>% 
  # Nimetään Long ja Lat X:ksi ja Y:ksi
  select(species, ID, year, X = Long, Y = Lat) %>% 
  # tidyr::spread toimii vain yhden sarakkeen kanssa, joten yhdistetään X ja Y
  # toistaiseksi samaan sarakkeeseen
  unite("XY", X, Y) %>%
  # Levitetään data yhdistelmäsarakkeen "XY" mukaan
  spread(year, XY) %>% 
  # Erotetaan 1960 X- ja Y-koordinaatit takaisin omiksi sarakkeikseen
  separate("60", c("X_60", "Y_60"), sep = "_") %>% 
  # Erotetaan 2010 X- ja Y-koordinaatit takaisin omiksi sarakkeikseen
  separate("10", c("X_10", "Y_10"), sep = "_") %>% 
  # Valitaan lopuksi sarakkeet oikeassa järjestyksessä
  select(species, ID, X_60, Y_60, X_10, Y_10) %>% 
  # Stringit numeroiksi
  mutate(X_60 = as.double(X_60), Y_60 = as.double(Y_60), 
         X_10 = as.double(X_10), Y_10 = as.double(Y_10))

# Grab the whole world, NOTE that you will also need package rworldextra
world.wgs84 <- getMap(resolution = "high")

world.laea <- spTransform(world.wgs84, crs.laea)

world.laea@data$border <- "darkgrey"

options(scipen = 9999)

plot(world.laea, xlim = c(4000000, 4700000), ylim = c(1500000, 5350000), 
     border = world.laea@data$border, axes = TRUE)

Arrows(x0 = sp_data_laea_change$X_60, y0 = sp_data_laea_change$Y_60,
       x1 = sp_data_laea_change$X_10, y1 = sp_data_laea_change$Y_10,
       arr.type = "triangle", arr.adj = 1, arr.length = 0.1)
text(x = jitter(rowMeans(cbind(sp_data_laea_change$X_60, sp_data_laea_change$X_10)), factor = 20),
     y = jitter(rowMeans(cbind(sp_data_laea_change$Y_60, sp_data_laea_change$Y_10)), factor = 20),
     labels = sp_data_laea_change$ID)