Skip to content

Instantly share code, notes, and snippets.

@mbacou

mbacou/snippets.R

Last active Nov 30, 2019
Embed
What would you like to do?
Working with attributes inside SpatialPolygonsDataFrame objects
## 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")
@sniffadog

This comment has been minimized.

Copy link

@sniffadog sniffadog commented Sep 20, 2015

Superb. Many thanks for this worked example. Very useful.

@cliffnyang

This comment has been minimized.

Copy link

@cliffnyang cliffnyang commented Mar 2, 2016

Very cool. Thanks for sharing!

@HanOostdijk

This comment has been minimized.

Copy link

@HanOostdijk HanOostdijk commented Mar 14, 2016

Thanks.

@Eflores89

This comment has been minimized.

Copy link

@Eflores89 Eflores89 commented Jan 31, 2017

Thanks! Quick question, how would you go about subsetting any one of the polygons inside the object? For example, eliminating ID number 1 (TZA)?

@vitaovtc123

This comment has been minimized.

Copy link

@vitaovtc123 vitaovtc123 commented Oct 12, 2017

I love you so mutch!!!!!!! Thank'ss

@realkrantz

This comment has been minimized.

Copy link

@realkrantz realkrantz commented Mar 1, 2018

This is very useful.

I want to assign a value from a variable " VAR1", contained in another dataset that also contains the variable " ID_2", to each value of the variable "ID_2" inside SpatialPolygonsDataFrame objects.

To accomplish this I replicated the code above except that I used inner_join with by = "ID_2" to add the new data (attribute) to the data of the SpatialPolygonsDataFrame object.

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.