Dealing with spatio-temporal data: sp::over function in R

Question

I have a dataset called carcass, which contains occurrences of livestock carcasses in Norway in the last 5 years (each raw corresponds to a single carcass, which has an ID number, along with longitude, latitude and year). I have also a shapefile called regions, a vector file over the different regions in Norway (every region has an ID number as well, specified in the column number_reg).

I managed to plot the map of the regions in R, and to plot the points corresponding to each carcass over the map.

regions<-readShapeSpatial("\\\\homer.uit.no/fma023/Desktop/Filippo Marolla/Datasets/rbd2012_siida.shp")
plot(regions)

carcass<-readShapePoints("\\\\homer.uit.no/fma023/Desktop/Filippo Marolla/Datasets/carcass_data_edit.shp")
plot(carcass, add=T, pch=20, cex=.5, col="red")

Then, I managed to count how many points fall within each district using the over function:

res<-table(over(carcass, regions)$number_reg

I also managed to plot the temporal trend of carcass occurrences for the entire country:

plot(table(carcass$YEAR), type="o", ylab="Carcass occurrences", main="Temporal trend of carcass occurrences in Norway")

Now I need to plot carcass occurrences over time for every single region, i.e. I want a graph with number of carcasses on y-axis and time (year) on x-axis for region 1, same graph for region 2, same graph for region 3, and so forth.

However, I am in trouble since the over function in R "(...) at the spatial location of object x retrieves the indexes or attributes from spatial object y", therefore in my res object I don't have the column year included (because year is in the carcass dataset, not in the regions dataset).

Answer 1

Here is a reproducible example using sp internal dataset meuse since I don't have access to your shapefiles:

r1 = cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409, 
             180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676, 
                                332618, 332413, 332349))
r2 = cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437, 
             179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683, 
                                        331133, 331623, 332152, 332357, 332373))
r3 = cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875, 
             179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
           c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004, 
             329783, 329665, 329720, 329933, 330478, 331062, 331086))
r4 = cbind(c(180304, 180403,179632,179420,180304),
           c(332791, 333204, 333635, 333058, 332791))

sr1=Polygons(list(Polygon(r1)),"r1")
sr2=Polygons(list(Polygon(r2)),"r2")
sr3=Polygons(list(Polygon(r3)),"r3")
sr4=Polygons(list(Polygon(r4)),"r4")
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf=SpatialPolygonsDataFrame(sr, data = data.frame(regions=1:4, row.names=c("r1","r2","r3","r4")))
plot(srdf)

# sample point dataset
data(meuse) 
#convert to spointdf
coordinates(meuse) = ~x+y
plot(meuse,add=T)

meuse_reg <- over(meuse,srdf[,1])
head(meuse_reg)
# add the landuse field
meuse_reg <- data.frame(meuse@data$landuse , meuse_reg)
# two-way contigency table for regions and landuse
plot(table(meuse_reg))

op <- par(mfrow=c(3,1))
plot(factor(meuse_reg[meuse_reg$region==1 ,1]))
plot(factor(meuse_reg[meuse_reg$region==2 ,1]))
plot(factor(meuse_reg[meuse_reg$region==3 ,1]))
par(op)

Check the code and adapt with your datasets. Here are the mappings from the above example to your data:

• yourData : myExample
• regions : regions
• carcass : meuse
• year : landuse (landuse and year are both categorical variables)
• number_reg : regs