2

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).

0

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
  • Happy to hear that. Please mark it as the correct answer as per tour – Farid Cheraghi May 27 '16 at 14:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.