WS 2016/2017

Prerequisites

setwd, library, data

setwd("YOUR-PATH")
raw_data <- read.table(file = "data/FS_UTM.csv",
                       header = TRUE,
                       sep = ";",
                       stringsAsFactors = FALSE,
                       quote = "")

4a. Check data structure

Do you remember the "4" issue from the previous session?! Let's check it again and solve the problem. I use a (s)apply way. You can chose the method you prefer.

First, let's check which columns have problems

#str(raw_data)
#sapply(raw_data[,4:20], function(x){x[x>1]})
apply(X=raw_data[,4:20], MARGIN = 2, FUN=function(x){x[x>1]})
## $settlement
## [1] 4
## 
## $open.settlement
## integer(0)
## 
## $cave.settlement
## integer(0)
## 
## $tell.settlement
## integer(0)
## 
## $mound
## integer(0)
## 
## $fortification
## integer(0)
## 
## $court
## integer(0)
## 
## $depot
## integer(0)
## 
## $cult.complex
## integer(0)
## 
## $mine
## integer(0)
## 
## $chruch
## integer(0)
## 
## $temple
## integer(0)
## 
## $villa.rustica
## integer(0)
## 
## $monastery
## integer(0)
## 
## $cemetery
## integer(0)
## 
## $necropolis
## integer(0)
## 
## $tumulus
## integer(0)

4a. Check data structure

Now, let's solve the problem and check again

raw_data$settlement[raw_data$settlement > 1] <- 1

sapply(raw_data[,4:20], function(x){x[x>1]})
## $settlement
## numeric(0)
## 
## $open.settlement
## integer(0)
## 
## $cave.settlement
## integer(0)
## 
## $tell.settlement
## integer(0)
## 
## $mound
## integer(0)
## 
## $fortification
## integer(0)
## 
## $court
## integer(0)
## 
## $depot
## integer(0)
## 
## $cult.complex
## integer(0)
## 
## $mine
## integer(0)
## 
## $chruch
## integer(0)
## 
## $temple
## integer(0)
## 
## $villa.rustica
## integer(0)
## 
## $monastery
## integer(0)
## 
## $cemetery
## integer(0)
## 
## $necropolis
## integer(0)
## 
## $tumulus
## integer(0)

4b. Check data structure

So, now it's time to have a look at the coordinates…

#str(raw_data)
raw_data[c(1:6),c(2,35,36)]
##                          Site.Name     xUTM    yUTM
## 1        Beba Veche - in the ditch 443157.0 5105906
## 2 Dudestii Vechi - Dragomirs' hill 459881.7 5095991
## 3               Nerau - Hunca Mare 462115.2 5088415
## 4             Nerau - Dogans' hill 466750.8 5089685
## 5          Cenad - Habitation Area 467929.9 5109462
## 6          Nerau - The south canal 468829.2 5092020

Question: Is there another way to create this overview?

4b. Check data structure

Yes.

#str(raw_data)
head(raw_data[,c(2,35,36)])
##                          Site.Name     xUTM    yUTM
## 1        Beba Veche - in the ditch 443157.0 5105906
## 2 Dudestii Vechi - Dragomirs' hill 459881.7 5095991
## 3               Nerau - Hunca Mare 462115.2 5088415
## 4             Nerau - Dogans' hill 466750.8 5089685
## 5          Cenad - Habitation Area 467929.9 5109462
## 6          Nerau - The south canal 468829.2 5092020

Create a spatial object

Some useful references

  • THE book that covers a broad range of topics:
  • Bivand, Roger S., Edzer J. Pebesma, and Virgilio Gómez-Rubio. 2008. Applied Spatial Data Analysis with R. New York: Springer.
  • Useful tutorial-like introduction in the handling and plotting of spatial data
  • Lovelace, Robin, and James Cheshire. 2015. “Introduction to Visualising Spatial Data in R.” link.

Spatial data libraries

Create spatial object

Create spatial data

sites <- SpatialPointsDataFrame(coords = cbind(raw_data$xUTM,raw_data$yUTM),
                                data = data.frame(raw_data[,2]),
                                proj4string = CRS("+init=epsg:32634")
                                )
class(sites)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

Exercise: Find other ways to create a spatial data object.

Create spatial data

sites <- raw_data
coordinates(sites) <- ~xUTM+yUTM
class(sites)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

Exercise: Check the structure of the sites object. What is missing?

is.projected(sites)
## [1] NA
proj4string(sites) <- CRS("+init=epsg:32634")
is.projected(sites)
## [1] TRUE

Plot it

Plot it 2.0

It's all about interaction, isn't it?

library(magrittr)
library(leaflet)

sites <- spTransform(sites, "+init=epsg:4326")

m1 <- leaflet() %>%
    addProviderTiles("Thunderforest.Landscape") %>%
        addMarkers(lng=sites@coords[,1],
                   lat=sites@coords[,2],
                   popup = sites@data$Site.Name
                   )
m1

Plot it 2.0

You can also change the map tile provider for a wide range of different maps. An overview can be found here: http://leaflet-extras.github.io/leaflet-providers/preview/index.html

library(magrittr)
library(leaflet)

sites <- spTransform(sites, "+init=epsg:4326")

m2 <- leaflet() %>%
    addProviderTiles("OpenStreetMap.HOT") %>%
        addMarkers(lng=sites@coords[,1],
                   lat=sites@coords[,2],
                   popup = sites@data$Site.Name
                   )
m2

Plot it 2.0

Another great package for interactive visualization: mapview

Question: What is the drawback of these interactive mapping approaches?

Export to Shapefile

writeOGR(obj = sites,
         dsn = "data",
         layer = "A_Sites",
         driver = "ESRI Shapefile",
         overwrite_layer = TRUE
         )
## Warning in writeOGR(obj = sites, dsn = "data", layer = "A_Sites", driver =
## "ESRI Shapefile", : Field names abbreviated for ESRI Shapefile driver
list.files("./data")
## [1] "A_Sites.dbf"    "A_Sites.prj"    "A_Sites.shp"    "A_Sites.shx"   
## [5] "FS_UTM.csv"     "FS_UTM.xlsx"    "srtm_41_03.tif"

Export to graphics file

not very useful; just to show how easy it is

png("data/A_Sites.png")
plot(sites)
dev.off()

## pdf, svg, postscript
?grDevices

?cairo

Play.