## Example procedure for working with data (attributes) inside SpatialPolygonsDataFrame objects | |
# Load my favorite libraries for that sort of work | |
library(data.table) | |
library(rgdal) | |
# Load 2 shapefiles that, say, we want to merge. | |
# Note that you can read pretty much any spatial format using readOGR(). | |
tza.l2 <- readOGR("./analysis/TZA-AC-07/maps", "tza-ac-07_L2") | |
tza.l1 <- readOGR("./analysis/TZA-AC-07/maps", "tza-ac-07_L1") | |
# Let's explore the 1st object (for information) | |
str(tza.l2) | |
# Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots | |
# ..@ data :'data.frame': 140 obs. of 12 variables: | |
# .. ..$ ISO : Factor w/ 1 level "TZA": 1 1 1 1 1 1 1 1 1 1 ... | |
# .. ..$ NAME_0 : Factor w/ 1 level "Tanzania": 1 1 1 1 1 1 1 1 1 1 ... | |
# .. ..$ ID_1 : int [1:140] 3052 3052 3052 3052 3053 3053 3053 3054 3054 3054 ... | |
# .. ..$ NAME_1 : Factor w/ 26 levels "Arusha","Dar-Es-Salaam",..: 1 1 1 1 2 2 2 3 3 3 ... | |
# .. ..$ ID_2 : int [1:140] 33658 33659 33661 33662 33664 33665 33666 33667 33668 33669 ... | |
# .. ..$ NAME_2 : Factor w/ 130 levels "Arumeru","Arusha",..: 1 2 76 96 20 40 123 13 14 44 ... | |
# .. ..$ Shape_Leng: num [1:140] 2.801 0.504 7.466 6.13 1.063 ... | |
# .. ..$ Shape_Area: num [1:140] 0.2321 0.00882 1.30787 1.25441 0.02747 ... | |
# .. ..$ svyL1Code : int [1:140] 2 2 2 2 7 7 7 1 1 1 ... | |
# .. ..$ svyL1Name : Factor w/ 26 levels "Arusha","Dar es salaam",..: 1 1 1 1 2 2 2 3 3 3 ... | |
# .. ..$ svyL2Code : int [1:140] 7 3 1 5 2 1 3 6 5 1 ... | |
# .. ..$ svyL2Name : Factor w/ 136 levels "Arusha Rural",..: 1 2 81 101 22 44 130 5 16 48 ... | |
# ..@ polygons :List of 140 | |
# .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots | |
# .. .. .. ..@ Polygons :List of 49 | |
# .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots | |
# .. .. .. .. .. .. ..@ labpt : num [1:2] 36.91 -2.99 | |
# .. .. .. .. .. .. ..@ area : num 1.63e-16 | |
# .. .. .. .. .. .. ..@ hole : logi FALSE | |
# .. .. .. .. .. .. ..@ ringDir: int 1 | |
# .. .. .. .. .. .. ..@ coords : num [1:5, 1:2] 36.9 36.9 36.9 36.9 36.9 ... | |
# .. .. .. .. ..[list output truncated] | |
# .. .. .. ..@ plotOrder: int [1:49] 31 48 23 45 40 38 15 37 16 20 ... | |
# .. .. .. ..@ labpt : num [1:2] 36.73 -3.43 | |
# .. .. .. ..@ ID : chr "0" | |
# .. .. .. ..@ area : num 0.12 | |
# .. .. [list output truncated] | |
# ..@ plotOrder : int [1:140] 92 47 110 60 114 73 97 13 116 100 ... | |
# ..@ bbox : num [1:2, 1:2] 29.327 -11.746 40.445 -0.982 | |
# .. ..- attr(*, "dimnames")=List of 2 | |
# .. .. ..$ : chr [1:2] "x" "y" | |
# .. .. ..$ : chr [1:2] "min" "max" | |
# ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots | |
# .. .. ..@ projargs: chr NA | |
# That object is made of 5 "slots" | |
slotNames(tza.l2) | |
# [1] "data" "polygons" "plotOrder" "bbox" "proj4string" | |
# We can access each slot (using "@") | |
class(tza.l2@data) | |
# [1] "data.frame" | |
class(tza.l2@polygons) | |
# [1] "list" | |
# Each polygon is itself a list with 5 slots (incl. all the vertice coordinates). | |
slotNames(tza.l2@polygons[[1]]) | |
# [1] "Polygons" "plotOrder" "labpt" "ID" "area" | |
# The ID slot is what ties the geometries and the attributes. | |
all(lapply(tza.l2@polygons, slot, "ID") == row.names(tza.l2)) | |
# [1] TRUE | |
# Note that like row.names, polygon IDs are character not integer. | |
summary(as.integer(lapply(tza.l2@polygons, slot, "ID"))) | |
# Min. 1st Qu. Median Mean 3rd Qu. Max. | |
# 0.0 34.8 69.5 69.5 104.0 139.0 | |
summary(as.integer(row.names(tza.l2)) | |
# Min. 1st Qu. Median Mean 3rd Qu. Max. | |
# 0.0 34.8 69.5 69.5 104.0 139.0 | |
# Save the data row.names into an explicit variable | |
tza.l2$rn <- row.names(tza.l2) | |
# Create temporary data tables to work on the attributes | |
tmp.l2 <- data.table(tza.l2@data) | |
tmp.l1 <- data.table(tza.l1@data) | |
# Say we want to merge data from table tmp.l1 into tmp.l2 | |
setkey(tmp.l2, svyL1Name, svyL2Code) | |
setkey(tmp.l1, region, district) | |
tmp.l2 <- tmp.l1[tmp.l2] | |
# Then let's re-attach the table to the original SpatialPolygonsDataFrame | |
# (preserving the original order of the row.names) | |
setkey(tmp.l2, rn) | |
tza.l2@data <- tmp.l2[row.names(tza)] | |
# Export back to shapefile (or to any spatial format) | |
writeOGR(tza.l2, "./analysis/TZA-AC-07/maps", "tza-ac-07_L2", driver="ESRI Shapefile") |
This comment has been minimized.
This comment has been minimized.
Very cool. Thanks for sharing! |
This comment has been minimized.
This comment has been minimized.
Thanks. |
This comment has been minimized.
This comment has been minimized.
Thanks! Quick question, how would you go about subsetting any one of the polygons inside the object? For example, eliminating ID number 1 (TZA)? |
This comment has been minimized.
This comment has been minimized.
I love you so mutch!!!!!!! Thank'ss |
This comment has been minimized.
This comment has been minimized.
This is very useful. I want to assign a value from a variable To accomplish this I replicated the code above except that I used But when I re-attach the table to the original SpatialPolygonsDataFrame I get Error in setkeyv (x, cols, verbose = verbose, physical = physical): x is not a data.table. I would appreciate any ideas. |
This comment has been minimized.
Superb. Many thanks for this worked example. Very useful.