Mapping Spatial Data

This markdown document shows some of the methods for importing and creating spatial dat in R and mapping these data in both static and interactive forms. There is a ton of great information in R mapping/GIS on the web. The authoritative source for keeping track of what exists is the CRAN Task View: Analysis of Spatial Data. Googling more specific packages will certainly get you some valuable information. Essentially, R is a fully functional GIS by any definition of the phrase. However, it is not a replacement for everything you may do in ArcGIS or QGIS, but it is a replacement for many tasks and an extension of these other software. Mapping in R really makes gains in spatial modeling, handling large data sets, repeatable tasks, and working spatial data into a larger analysis

The example here uses the Michelsberg data set from the archdata package. These data include the location and attributes of vessels gathered from 108 Michelsberg period features. The code below walks through how to make these data into spatial objects, how to import a shapefile, create a static map, conduct spatial analysis on a raster, and finally map these layers in an interactive map.

Data details:

A sites by types table of abundance data on vessel types in archaeological features of the Younger Neolithic Michelsberg Culture from Belgium, France and Germany by Birgit Höhn (2002).

Load libraries

library("archdata") # where the data comes from
library("sp") # for spatial object
library("leaflet") # for interactive mapping
library("rgdal") # for geographic transformations
library("rgeos") # for geographic operations
library("spatstat")  # Needed for the dirichlet tesselation function
library("rworldmap") # library rworldmap provides different types of global maps, e.g:
library("ggmap") # for a static map example
library("gstat") # for IDW model
library("leaflet") # for interactive mapping
library("mapview") # for interactive mapping test with mapview(breweries91)
library("raster") # for raster operations

Invoke data and note structure

### bring in data from archdata package
data(Michelsberg)
# give a shorter name, less typing...
mb <- Michelsberg
# take a look at structure of data
str(mb)
'data.frame':   109 obs. of  42 variables:
 $ id          : int  2 1 3 4 5 6 7 8 9 10 ...
 $ site_name   : chr  "achenheim" "achenheim" "didenheim" "entzheim" ...
 $ catalogue_nr: int  1 1 2 3 4 5 6 8 8 9 ...
 $ feature_nr  : num  1.1 1.2 2 3 4 5 6 8.1 8.2 9 ...
 $ to3         : int  0 0 0 0 0 0 0 1 0 0 ...
 $ f4          : int  0 0 0 0 0 0 0 0 0 0 ...
 $ b2          : int  0 0 1 0 3 0 0 0 0 0 ...
 $ to2         : int  0 0 0 0 0 0 0 0 0 0 ...
 $ b3          : int  0 0 0 0 0 0 0 0 0 0 ...
 $ b7          : int  0 0 0 0 1 0 0 0 0 0 ...
 $ kw5         : int  0 0 1 0 0 0 0 0 2 0 ...
 $ vg1         : int  0 0 1 0 8 0 0 2 1 0 ...
 $ vg2         : int  0 0 0 0 0 0 0 0 0 0 ...
 $ t4a         : int  0 0 0 0 0 0 0 0 0 0 ...
 $ kw2         : int  0 1 1 0 1 2 1 1 5 0 ...
 $ kw4         : int  0 0 0 0 1 0 0 1 1 0 ...
 $ b5          : int  0 0 2 0 0 0 0 0 0 0 ...
 $ t3b         : int  0 0 0 0 0 0 0 0 0 0 ...
 $ f3          : int  0 0 0 0 0 0 0 0 0 0 ...
 $ kw3         : int  0 2 1 0 0 2 1 0 4 0 ...
 $ kw1         : int  0 0 0 0 0 0 0 0 0 0 ...
 $ b6          : int  0 0 1 0 0 0 0 0 0 0 ...
 $ to1         : int  0 0 0 0 0 0 0 0 0 0 ...
 $ b1          : int  0 1 0 0 0 0 1 0 0 0 ...
 $ t3a         : int  0 0 0 0 0 0 1 0 0 0 ...
 $ vg4         : int  0 0 2 0 1 0 0 0 0 1 ...
 $ ks2         : num  0 2 0 0 0 0 0.5 0 0 0 ...
 $ ks1         : num  0 0 0 0 0 0 0.5 0 0 0 ...
 $ t2b         : int  0 1 0 0 0 3 2 0 0 0 ...
 $ f2          : int  0 0 0 0 0 0 0 0 0 0 ...
 $ bs3         : int  0 0 0 0 0 0 0 0 0 0 ...
 $ t2a         : int  0 1 0 0 0 3 1 0 0 1 ...
 $ bs2         : int  0 0 0 0 0 0 0 0 0 0 ...
 $ b4          : int  1 0 0 1 0 0 0 0 0 0 ...
 $ bs1         : int  1 0 0 1 0 0 0 0 0 0 ...
 $ f1          : int  1 0 0 0 0 0 0 0 0 1 ...
 $ t1b         : int  0 0 0 2 0 1 0 0 0 1 ...
 $ vg3         : int  2 0 0 0 0 0 0 0 0 2 ...
 $ t1a         : num  0 0 0 2 0 0 0 0 0 0 ...
 $ mbk_phase   : Ord.factor w/ 11 levels "I"<"I/II"<"II"<..: 2 5 10 1 10 5 5 10 10 3 ...
 $ x_utm32n    : num  398587 398587 372268 398856 399555 ...
 $ y_utm32n    : num  5381681 5381681 5286029 5376116 5373880 ...
 - attr(*, "typological_key")= chr  "b_ = beaker (Becher)" "bs_ = globular bowl (beckenförmige Schüssel)" "f_ = bottle (Flasche)" "ks_ = conical bowl (Konische Schüssel)" ...
 - attr(*, "references")=List of 2
  ..$ reference_1_to_108: chr "Höhn, B. 2002. Die Michelsberger Kultur in der Wetterau. Universitätsforschungen zur Prähistorischen Archäologie. Bonn: Habelt."| __truncated__
  ..$ reference_109     : chr "Mischka, D., Roth, G. and K. Struckmeyer 2015. Michelsberg and Oxie in contact next to the Baltic Sea. In: Kr. Brink et al. (ed"| __truncated__
 - attr(*, "projection")= chr "UTM32n_WGS84"
# take a look at type and range of variables
summary(mb)
       id       site_name          catalogue_nr     feature_nr   
 Min.   :  1   Length:109         Min.   : 1.00   Min.   : 1.10  
 1st Qu.: 28   Class :character   1st Qu.:19.00   1st Qu.:19.38  
 Median : 55   Mode  :character   Median :28.00   Median :28.15  
 Mean   : 55                      Mean   :33.94   Mean   :34.09  
 3rd Qu.: 82                      3rd Qu.:51.25   3rd Qu.:51.25  
 Max.   :109                      Max.   :69.00   Max.   :69.00  
                                  NA's   :1       NA's   :1      
      to3              f4               b2               to2        
 Min.   :0.000   Min.   :0.0000   Min.   : 0.0000   Min.   :0.0000  
 1st Qu.:0.000   1st Qu.:0.0000   1st Qu.: 0.0000   1st Qu.:0.0000  
 Median :0.000   Median :0.0000   Median : 0.0000   Median :0.0000  
 Mean   :0.156   Mean   :0.1651   Mean   : 0.6789   Mean   :0.0367  
 3rd Qu.:0.000   3rd Qu.:0.0000   3rd Qu.: 1.0000   3rd Qu.:0.0000  
 Max.   :6.000   Max.   :3.0000   Max.   :18.0000   Max.   :2.0000  
                                                                    
       b3                b7              kw5               vg1         
 Min.   : 0.0000   Min.   :0.0000   Min.   : 0.0000   Min.   : 0.0000  
 1st Qu.: 0.0000   1st Qu.:0.0000   1st Qu.: 0.0000   1st Qu.: 0.0000  
 Median : 0.0000   Median :0.0000   Median : 0.0000   Median : 0.0000  
 Mean   : 0.3303   Mean   :0.2018   Mean   : 0.6789   Mean   : 0.5596  
 3rd Qu.: 0.0000   3rd Qu.:0.0000   3rd Qu.: 1.0000   3rd Qu.: 1.0000  
 Max.   :15.0000   Max.   :4.0000   Max.   :22.0000   Max.   :11.0000  
                                                                       
      vg2              t4a              kw2               kw4        
 Min.   :0.0000   Min.   :0.0000   Min.   : 0.0000   Min.   : 0.000  
 1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.: 0.0000   1st Qu.: 0.000  
 Median :0.0000   Median :0.0000   Median : 0.0000   Median : 0.000  
 Mean   :0.2202   Mean   :0.1743   Mean   : 0.8991   Mean   : 0.945  
 3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.: 1.0000   3rd Qu.: 0.000  
 Max.   :4.0000   Max.   :4.0000   Max.   :25.0000   Max.   :24.000  
                                                                     
       b5              t3b               f3               kw3        
 Min.   :0.0000   Min.   :0.0000   Min.   :0.00000   Min.   : 0.000  
 1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.: 0.000  
 Median :0.0000   Median :0.0000   Median :0.00000   Median : 1.000  
 Mean   :0.1743   Mean   :0.2661   Mean   :0.08257   Mean   : 2.229  
 3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.: 2.000  
 Max.   :2.0000   Max.   :4.0000   Max.   :2.00000   Max.   :41.000  
                                                                     
      kw1               b6               to1                b1        
 Min.   : 0.000   Min.   :0.00000   Min.   :0.00000   Min.   :0.0000  
 1st Qu.: 0.000   1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.0000  
 Median : 0.000   Median :0.00000   Median :0.00000   Median :0.0000  
 Mean   : 1.092   Mean   :0.08257   Mean   :0.05505   Mean   :0.2385  
 3rd Qu.: 1.000   3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.0000  
 Max.   :18.000   Max.   :1.00000   Max.   :2.00000   Max.   :3.0000  
                                                                      
      t3a              vg4              ks2              ks1        
 Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
 1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
 Median :0.0000   Median :0.0000   Median :0.0000   Median :0.0000  
 Mean   :0.2477   Mean   :0.3211   Mean   :0.2248   Mean   :0.1605  
 3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000  
 Max.   :4.0000   Max.   :4.0000   Max.   :4.5000   Max.   :5.5000  
                                                                    
      t2b               f2              bs3              t2a        
 Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
 1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
 Median :0.0000   Median :0.0000   Median :0.0000   Median :0.0000  
 Mean   :0.5505   Mean   :0.1927   Mean   :0.2661   Mean   :0.6514  
 3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:1.0000  
 Max.   :6.0000   Max.   :6.0000   Max.   :8.0000   Max.   :6.0000  
                                                                    
      bs2                b4              bs1               f1        
 Min.   : 0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
 1st Qu.: 0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
 Median : 0.0000   Median :0.0000   Median :0.0000   Median :0.0000  
 Mean   : 0.7982   Mean   :0.1284   Mean   :0.3578   Mean   :0.2018  
 3rd Qu.: 0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000  
 Max.   :18.0000   Max.   :4.0000   Max.   :8.0000   Max.   :3.0000  
                                                                     
      t1b               vg3              t1a           mbk_phase 
 Min.   : 0.0000   Min.   :0.0000   Min.   :0.0000   III    :24  
 1st Qu.: 0.0000   1st Qu.:0.0000   1st Qu.:0.0000   V      :16  
 Median : 0.0000   Median :0.0000   Median :0.0000   II     :14  
 Mean   : 0.7248   Mean   :0.3303   Mean   :0.3394   Munz   :14  
 3rd Qu.: 1.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   II/III : 9  
 Max.   :12.0000   Max.   :5.0000   Max.   :8.0000   IV     : 9  
                                                     (Other):23  
    x_utm32n         y_utm32n      
 Min.   :116831   Min.   :5282624  
 1st Qu.:392708   1st Qu.:5408467  
 Median :467581   Median :5482002  
 Mean   :436610   Mean   :5493652  
 3rd Qu.:511310   3rd Qu.:5576821  
 Max.   :627228   Max.   :6010300

Making data spatial

create objects for coordinate referrence systems

## Proj4string from http://spatialreference.org/
# UTM 1983 zone 32 North - Michelsberg data
UTM32n <- CRS("+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
# World Geographic System 1984 (lat/long) - mapping
WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")

Create a spatialpointsdataframe object from our data

## pass long, lat, and data to function.  
## proj4string = the proper CRS as described in `archdata` package reference
mb_spat <-  SpatialPointsDataFrame(coords = mb[,c("x_utm32n", "y_utm32n")], 
                                   data = mb, 
                                   proj4string = UTM32n)

Static mapping with base plot()

First make a static map with the base plot() function

## most basic plot (points)
plot(mb_spat)

Putting it in geographic context

Bring in country boundary data from the rworldmap package, reproject the mb data to WGS84 and use base plot() again to plot vessel locations, country boundaries, and labels.

# get country boundary data
data(countriesLow)
# transform Michelsberg site data from UTM83n to WGS84
mb_spat_WGS84 <- spTransform(mb_spat, WGS84)
# plot the data and boundaries, then lable
plot(mb_spat_WGS84, pch = 20, col = "red")
# 'add' tell plot to add to the previous plot
plot(countriesLow, add = TRUE)
labelCountries()

ggmap example: Static map with nice basemap

The next approach is to use the ggmap package to create nice base maps and use the ggplot framework for mapping.

# grab WGS84 coordinates
WGS84_coords <- coordinates(mb_spat_WGS84)
# append WGS84 coords onto original data
mb$x_WGS84 <- WGS84_coords[,1]
mb$y_WGS84 <- WGS84_coords[,2]
# get a base map centered on Frankfurt, Germany at a certain zoom
base_map <- get_map(location = "Frankfurt, Germany", zoom = 6, color = "bw")
Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=Frankfurt,+Germany&zoom=6&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Frankfurt,%20Germany&sensor=false
# use ggplot code to build out layers of map:
# basemap, add points, change color scale, set plot theme
ggmap(base_map, extent = "normal") +
  geom_point(aes(x = x_WGS84, y = y_WGS84, color = t2a), data = mb, alpha = .8, size = 4) +
  scale_color_distiller(palette = "PuOr", direction = 1) +
  theme_bw(16)

Loading and adding a shapefile to the plot

Loading a shapefile is pretty easy with the readOGR() function of the rgdal package. If the shapefile has a crs assigned to it, readOGR() will detect it and assign it to the spatial object in R; very handy. For this demo, we load a shapefile for the border of Germany [in the /data/ folder on GitHub]. Next we use gSimplify() to reduce the number of edges to speed up the rendering. Finally, base plot() is used to take a look.

wd <- "/Users/mattharris/Documents/R_Local/SAA_R_Demo_2016/SAA_R_Demo_2016"
germany <- readOGR(paste0(wd, "/data/DEU_adm0.shp"), layer = "DEU_adm0")
OGR data source with driver: ESRI Shapefile 
Source: "/Users/mattharris/Documents/R_Local/SAA_R_Demo_2016/SAA_R_Demo_2016/data/DEU_adm0.shp", layer: "DEU_adm0"
with 1 features
It has 70 fields
germany <- gSimplify(germany, tol=0.001, topologyPreserve=TRUE)
plot(germany)

Fortifying for ggplot

To use the imported shapefile in ggmap, we need the fortify() function of the ggplot2 package. The fortify() function converts the spatial data into a dataframe suitable for plotting in ggmap or ggplot. As with a typical GIS, the plot is build by layers, first the geom_polygon() of the boundary, then the geom_point() of the vessel locations, and finally some color and display settings.

Note:

the ggplot2 function of coord_map() is used here to “zoom” the map into the extents of the data. There are other functions such as ylim() and scale_y_continuous() that control plot extents, but you will need coord_map() to do so with map data. Otherwise, data is cut off or shifted out of alignment.

germany_fort <- ggplot2::fortify(germany)
# basemap with points and SHP
ggmap(base_map, extent = "normal") +
  geom_polygon(data = germany, aes(x=long, y=lat, group=group), 
               fill = "transparent", color = "gray10", size = 1) +
  geom_point(aes(x = x_WGS84, y = y_WGS84, color = t2a), 
             data = mb, alpha = .8, size = 4) +
  scale_color_distiller(palette = "PuOr", direction = 1) +
  coord_map(ylim = c(min(mb$y_WGS84), max(mb$y_WGS84)),
              xlim = c(min(mb$x_WGS84), max(mb$x_WGS84))) +
  theme_bw(16)

Spatial Analysis

The first foray into spatial analysis happens directly in out calls to ggmap(). The ggplot2 framework is so extensive, it even has a number of statistical functions that can manipulate the data directly in the plotting function. In this case, stat_bin2d() is used to aggregate the number of vessels into aerial units.

Binned density

## use ggplot functions to creat site density grid
ggmap(base_map, extent = "normal") +
  stat_bin2d(aes(x = x_WGS84, y = y_WGS84, color = t2a),
  size = 1, bins = 20, alpha = 0.9, data = mb)

Smooth density

Using the stat_density2d() function

ggmap(base_map, extent = "normal") +
  stat_density2d(aes(x = x_WGS84, y = y_WGS84, fill = ..level.., alpha = ..level..),
  size = 1, bins = 20, data = mb,
  geom = "polygon")

Faceted density

Finally, we create a multi-panel, or faceted, map that shows the data distribution of data filtered on the number of occurrences of vessel type "t2a" #### Note: this can take a while to render!

mb$group <- ifelse(mb$t2a <= 3, "t2a <= 3", "t2a > 3")
ggmap(base_map, extent = "normal") +
  stat_density2d(aes(x = x_WGS84, y = y_WGS84, fill = ..level.., alpha = ..level..),
                 size = 1, bins = 20, data = mb,
                 geom = "polygon") +
  scale_fill_distiller(palette = "PuOr", direction = 1) +
  theme_bw() +
  facet_wrap(~ group)

Analysis with gstat

Next we perform a bit more in depth spatial analysis outside of the ggmap() call and then layer it back in.

This code uses the gstat package to perform an Inverse Distance Weighting (IDW) interpolation on our data (sort of… a new vessel type is created and given a fake distribution that looks better for demonstration purposes).

Build a grid to the geographic limits of the data

# get the min/max range for lat/long to make an empty grid 
x.range <- as.numeric(c(min(mb$x_WGS84), max(mb$x_WGS84)))  # min/max longitude of the interpolation area
y.range <- as.numeric(c(min(mb$y_WGS84), max(mb$y_WGS84)))  # min/max latitude of the interpolation area  
# from the range, exapnd the coordinates to make a regular grid
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.075), 
                   y = seq(from = y.range[1], to = y.range[2], by = 0.075))  # expand

make the grid into spatial coords, grid, assign WGS84 CRS

coordinates(grd) <- ~x + y
sp::gridded(grd) <- TRUE
proj4string(grd) <- WGS84
# quick plot to see what we have
plot(grd, cex = 1.5, col = "grey")
points(mb_spat_WGS84, pch = 1, col = "red", cex = 1)

simulate a new vesse type with a nicely defined large scale geographic distribution trend; call it new_type

### Because most attributes in example data don't have a great distribution
### I am simulating a made up attribute so tha mapping looks nicer
## simulate 'new_type' with spatial trend
# finds mean of site locations,
# assigns a value that increases from lower-left to upper-right
x_mean <- mean(mb$x_WGS84)
y_mean <- mean(mb$y_WGS84)
x_dist <- mb$x_WGS84 - x_mean
y_dist <- mb$y_WGS84 - y_mean
mb$xy_dist <- x_dist + y_dist
mb$pre_sort_id <- seq(1,nrow(mb))
mb <- mb[order(mb$xy_dist),]
mb$new_type <- seq(1,nrow(mb))*2
mb <- mb[order(mb$pre_sort_id),]
mb_spat_WGS84$new_type <- mb$new_type

Perform IDW interpolation

The gstat::IDW() function performs the interpolation and then the results are converted into a data.frame using as.data.frame().

### IDW interpolation of 'new_type' using gstat package
idw <- gstat::idw(formula = new_type ~ 1, locations = mb_spat_WGS84, newdata = grd)  # apply idw model for the data
[inverse distance weighted interpolation]
# grab output of IDW for plotting
idw.output = as.data.frame(idw)  # output is defined as a data table
# set the names of the idw.output columns
# basic ggplot using geom_tile to display our interpolated grid within no map
ggplot() + 
  geom_tile(data = idw.output, aes(x = x, y = y, fill = var1.pred)) + 
  geom_point(data = mb, aes(x = x_WGS84, y = y_WGS84), shape = 21, color = "red") +
  scale_fill_distiller(palette = "PuOr", direction = 1) +
  theme_bw()

Add IDW results to ggmap()

To add the IDW results to the ggmap, the geom_tile() function is used to create a raster display of the interpolated values, the result of the ggmap() call is about the same.

# use ggmap as earlier to display interpolation in geographic space
# Can take a LONG time to render (~1 minute)
ggmap(base_map, extent = "normal") +
  geom_tile(data = idw.output, aes(x = x, y = y, fill = var1.pred), alpha = 0.75) + 
  geom_point(data = mb, aes(x = x_WGS84, y = y_WGS84), shape = 19, color = "red") +
  scale_fill_distiller(palette = "PuOr", direction = 1) +
  theme_bw()

Interactive Mapping

Static maps are fantastic, but sometimes “slippy” interactive maps are great for presenting and exploring data. The following code uses the mapview package to create the interactive maps. This is a rapidly developing area for R mapping and other solutions currently exist, including leaflet.

The process of making an interactive mapview map is pretty darn easy as we will see. The first example is created simply by calling the mapview() function on the germany spatial object (imported from a shapefile)

mapview(germany)

Multi-layer map

Adding layers is a sample as creating a call to the mapview() function for each layer, and then adding them together. Styling for each layer is done within the call to mapview().

maplayer1 <- mapview(germany, 
                     color = "gray10", 
                     alpha.regions = 0)
maplayer2 <- mapview(mb_spat_WGS84)
maplayer1 + maplayer2

Adding in the IDW raster

The first step to incorporate the IDW results into the interactive map is to turn it into a spatial object. This is done with the function rasterFromXYZ() from the raster package. Plot it using ssplot() to make sure it looks reasonable.

## use raster packaage to turn IDW into raster object
idw_raster <- rasterFromXYZ(idw.output[,1:3], crs = WGS84)
# plot to see what happens
spplot(idw_raster)

Clip IDW raster

For this analysis, we only want to show the interpolation results within the area of which we have data points. To do so, we will simply clip the raster. To do so: 1. Create clipping polygon from UTM projected feature location data + Find the envelope of the data points with gConvexHull() + Buffer that envelope to make our study region with gBuffer() + reproject to WGS84 for mapping 2. Clip IDW raster + use mask() function in raster package to clip raster with buffer 3. Create contours + use rasterToContour() function to create contours of IDW results

## To clip raster to data region, creat convex hull, buffer, and mask (clip) raster
## perform convex hull and buffer operations on projected UTM32n data
mb_hull <- gConvexHull(mb_spat) # convex hull
mb_buff <- gBuffer(mb_hull, width = 25000) # arbitrary 25km buffer
# tranforms back to WGS84
mb_buff_WGS84 <- spTransform(mb_buff, WGS84)
# plot to see results
plot(mb_buff_WGS84)
points(mb_spat_WGS84)

# mask IDW raster to confine to region that we have sites for
idw_raster_crop <- mask(idw_raster, mb_buff_WGS84)
# create contours of the IDW
idw_contour <- rasterToContour(idw_raster_crop, nlevels = 15)
# plot to see results
spplot(idw_raster_crop)

Final map

Finally, all of the layers (Germany border, vessel locations, IDW raster, and contours) add added together as mapview objects.

### Use mapview package to create layered interactive map with
### basemap, data points, interpolated raster, and contours
## made in three seperate layers
m1 <- mapview(germany, color = "black", alpha.regions = 0)
m2 <- mapview(mb_spat_WGS84, legend = TRUE, zcol = "new_type")
m3 <- mapview(idw_raster_crop, alpha.regions = 0.50, na.color = "transparent")
m4 <- mapview(idw_contour, lwd = 1.5, color = "gray40")
# combine layers to plot
m1 + m2 + m3 + m4

Data Details

Description

A sites by types table of abundance data on vessel types in archaeological features of the Younger Neolithic Michelsberg Culture from Belgium, France and Germany by Birgit Höhn (2002).

Details

Höhn (2002) recorded pottery vessel shapes from 108 archaeological features (pits, ditches etc.) from 69 sites of the Central European Younger Neolithic Michelsberg Culture (MBK; 4350-3500 BC) following Lüning’s (1967) typology. Her correspondence analysis of the abundance data (columns 5 to 39) exhibits a pronounced Guttman effect or arch, suggesting the data set is structured by a time gradient. Recently Mischka et al. (2015) projected an 109th Michelsberg assemblage, Flintbek LA48, a pit with Michelsberg pottery from a North German site of the Funnel Beaker Cul- ture (TRB), as a supplementary row into the existing chronology thereby connecting the relative chronologies of TRB and MBK. The data frame contains as attributes the references for the data, a typological key and the map projection. Note that ambiguous fragments of conical bowls (ks1 and ks2) are assigned as 0.5 to each of the two types resulting also in positive entries suitable to analysis by CA.

Source

Höhn, B. 2002. Die Michelsberger Kultur in der Wetterau. Universitätsforschungen zur prähis- torischen Archäologie 87. Bonn: Habelt. Mischka, D., Roth, G. and K. Struckmeyer 2015. Michelsberg and Oxie in contact next to the Baltic Sea. In: Neolithic Diversities. Perspectives from a conference in Lund, Sweden. Acta Archaeologica Lundensia Ser. 8, No. 65, edited by Kr. Brink et al., pp 241–250. Lüning, J. 1967. Die Michelsberger Kultur: Ihre Funde in zeitlicher und räumlicher Gliederung. Berichte der Römisch-Germanischen Kommission 48, 1-350.