Simple_site_grid
Matthew Harris
April 19, 2017
Example
Simple example of simulated archaeological excavation data and the spatialPixels and spatialPolygons dataframe classes.
libraries
library("dplyr") # data manipulation
library("tidyr") # data manipulation
library("sp") # spatial classes
library("mapview") # interactive mappingCoordinate System
Using something local to me for no real reason…
UTM18N83 <- "+proj=utm +zone=18 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"Site Grid
We will use the full site grid even though ~60% does not have data. This will help us frame the data and make the sp classes. I made a square 19x19 grid here to start just to simplify the simulation of your data. In the first bit, to site_grid I assign the row/col letters/numbers and also their corresponding UTM coords
In the second bit, I create grid_data by nesting and expanding the coordinates into all combinations of them. You can arrive at your data in this structure a number of way; this just made it simple to simulate. In that call, I also concatenate fields into the unit names and coordinate pairs (probably use that later somewhere). Finally, I filter out where the y coord is greater than 17 so I arrive at the same 19x17 grid you described
site_grid <- data.frame( x = letters[1:19],
y = seq(1,19,1),
x_coord = seq(481744,481762,1),
y_coord = seq(4425933,4425951,1))
grid_data <- site_grid %>%
tidyr::expand(nesting(x,x_coord),nesting(y,y_coord)) %>%
mutate(name = paste(x,y, sep = ''),
coord = paste(x_coord, y_coord, sep = ',')) %>%
filter(y <= 17)
head(grid_data)## # A tibble: 6 × 6
## x x_coord y y_coord name coord
## <fctr> <dbl> <dbl> <dbl> <chr> <chr>
## 1 a 481744 1 4425933 a1 481744,4425933
## 2 a 481744 2 4425934 a2 481744,4425934
## 3 a 481744 3 4425935 a3 481744,4425935
## 4 a 481744 4 4425936 a4 481744,4425936
## 5 a 481744 5 4425937 a5 481744,4425937
## 6 a 481744 6 4425938 a6 481744,4425938Field Data
To simulate excavated field data, I take the 19x17 site grid, simulate a count of fictitious artifacts for each cell and then select 40% of the cells at random. This equates to having excavated 40% of the site.
field_data <- grid_data %>%
mutate(art_cnt = floor(rnorm(nrow(grid_data),20,5))) %>%
sample_frac(0.40) %>%
dplyr::select(name, art_cnt)
head(field_data)## # A tibble: 6 × 2
## name art_cnt
## <chr> <dbl>
## 1 n11 25
## 2 n2 19
## 3 r5 19
## 4 h15 29
## 5 g14 23
## 6 p7 15Plotting Data
In order to plot or work with the field_data as an SP object, the artifacts counts need to be joined to the full site grid where a unit name matches a grid cell name (e.g. “L12”). The resulting join returns all of the grid cells and a column for artifact count (art_cnt) that is a number where an excavated unit matches the grid (they all hopefully do match) and an NA where the cell is unexcavated.
plot_data <- left_join(grid_data, field_data, by = "name")
head(plot_data)## # A tibble: 6 × 7
## x x_coord y y_coord name coord art_cnt
## <fctr> <dbl> <dbl> <dbl> <chr> <chr> <dbl>
## 1 a 481744 1 4425933 a1 481744,4425933 NA
## 2 a 481744 2 4425934 a2 481744,4425934 NA
## 3 a 481744 3 4425935 a3 481744,4425935 NA
## 4 a 481744 4 4425936 a4 481744,4425936 20
## 5 a 481744 5 4425937 a5 481744,4425937 NA
## 6 a 481744 6 4425938 a6 481744,4425938 36spatialPixelsDataFrame
So your question was specific to spatialP*DataFrames, but of course there are a number of ways to visualize and analyze these data; I’ll just look at getting it into the SP class for now.
This is as simple as calling the SpatialPixelsDataFrame() function and using the spatial coordinates, data, and projection as the parameters.
SPixels <- SpatialPixelsDataFrame(points = plot_data[,c("x_coord","y_coord")],
data = plot_data[,"art_cnt"],
proj4string = CRS(UTM18N83))
spplot(SPixels)
spatialPolygonsDataFrame
This can be made into a SpatialPolygonsDataFrame with a simple conversion.
SPolygons = as(SPixels, "SpatialPolygonsDataFrame") ## Warning: Setting row names on a tibble is deprecated.spplot(SPolygons) 
mapview Example
And to show that this is a proper spatial object, mapview() is used here to plot it at the arbitrary coordinates I had give the grid. Note: as this is only a 19x17-meter grid, the zoom level of this map makes it fairly useless. However, if you zoom out a bit, you will see the grid and the background are in real geographic space. You can expand on this a great deal for more interactive visualization and inspection.
mapview(SPixels) sessionInfo()
sessionInfo()## R version 3.3.2 (2016-10-31)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 7 x64 (build 7601) Service Pack 1
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] mapview_1.2.0 leaflet_1.1.0 sp_1.2-4 tidyr_0.6.1 dplyr_0.5.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.10 plyr_1.8.4 RColorBrewer_1.1-2
## [4] base64enc_0.1-3 R.methodsS3_1.7.1 R.utils_2.5.0
## [7] iterators_1.0.8 tools_3.3.2 gdalUtils_2.0.1.7
## [10] digest_0.6.12 jsonlite_1.2 viridisLite_0.1.3
## [13] satellite_0.2.0 evaluate_0.10 tibble_1.2
## [16] lattice_0.20-34 png_0.1-7 foreach_1.4.3
## [19] shiny_1.0.0 DBI_0.5-1 crosstalk_1.0.0
## [22] yaml_2.1.14 rgdal_1.2-5 stringr_1.2.0
## [25] knitr_1.15.1 raster_2.5-8 htmlwidgets_0.8
## [28] webshot_0.4.0 stats4_3.3.2 rprojroot_1.2
## [31] grid_3.3.2 R6_2.2.0 rmarkdown_1.3
## [34] latticeExtra_0.6-28 magrittr_1.5 scales_0.4.1
## [37] backports_1.0.5 codetools_0.2-15 htmltools_0.3.5
## [40] assertthat_0.1 colorspace_1.3-2 mime_0.5
## [43] xtable_1.8-2 httpuv_1.3.3 stringi_1.1.2
## [46] lazyeval_0.2.0 munsell_0.4.3 R.oo_1.21.0