WS 2016/2017
setwd("YOUR-PATH")
raw_data <- read.table(file = "data/FS_UTM.csv", header = TRUE, sep = ";", stringsAsFactors = FALSE, quote = "")
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)
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)
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?
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
library(sp) library(rgdal)
As always there are numerous ways to accomplish this. See https://cran.r-project.org/web/packages/sp/vignettes/intro_sp.pdf
There are some things we need to know in advance:
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.
sites <- raw_data coordinates(sites) <- ~xUTM+yUTM class(sites)
## [1] "SpatialPointsDataFrame" ## attr(,"package") ## [1] "sp"
Exercise: Check the str
ucture of the sites
object. What is missing?
is.projected(sites)
## [1] NA
proj4string(sites) <- CRS("+init=epsg:32634") is.projected(sites)
## [1] TRUE
plot(sites)
Find further ways to plot the points.
A collection of examples: https://pakillo.github.io/R-GIS-tutorial/
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
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
Another great package for interactive visualization: mapview
Question: What is the drawback of these interactive mapping approaches?
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"
not very useful; just to show how easy it is
png("data/A_Sites.png") plot(sites) dev.off() ## pdf, svg, postscript ?grDevices ?cairo