0

I have a spatial point data frame with two columns with attributes, one for months, another for species (and off-course the lat, long for those points).

long=c(-1.747429,-1.519554,0.430455,-2.049941,-1.288655,1.311327,-2.217094,0.718583,-1.646491)

lat= c(54.407632,53.369744,51.519982,52.347591,54.530768,51.12631,53.430349,51.552895,55.028696)

species = c('d','c','d','d','d','c','c','c','c')

month = c(1,2,3,2,3,4,3,2,3)

I also have a spatial polygon of my country projected in the same CRS as my points.

I want to create a loop that will give me the count of points per month and species and that output to be in separate data frames.

But so far I only know how to do count points in polygon overall

Df2<-over(Df, poly)

df2Freq<-count(Df2, "NAME_1")

What I want is a loop that will do by month and species and produce multiple data frames with outputs

1

I will provide you a reproducible example. First, I'll create a point object and a polygon. The point object has month and species field:

library(sp)
library(raster)

point <- new("SpatialPointsDataFrame"
             , data = structure(list(month = c(1, 1, 1, 1, 2, 2, 2, 3, 3, 3), species = c(1, 
                                                                                          2, 3, 1, 2, 3, 3, 1, 3, 2)), .Names = c("month", "species"), row.names = c(NA, 
                                                                                                                                                                     -10L), class = "data.frame")
             , coords.nrs = numeric(0)
             , coords = structure(c(-71.2611979005976, -71.2613792984669, -71.262302778529, 
                                    -71.2594333940505, -71.2594828661966, -71.2629624071448, -71.263556072899, 
                                    -71.2606866884204, -71.2589881447348, -71.258658330427, -29.9447590462946, 
                                    -29.9441159083943, -29.9436541683632, -29.9428131418781, -29.9417577360929, 
                                    -29.945484637772, -29.9446271205715, -29.9421040411162, -29.941972115393, 
                                    -29.9423019297009), .Dim = c(10L, 2L), .Dimnames = list(NULL, 
                                                                                            c("coords.x1", "coords.x2")))
             , bbox = structure(c(-71.263556072899, -29.945484637772, -71.258658330427, 
                                  -29.9417577360929), .Dim = c(2L, 2L), .Dimnames = list(c("coords.x1", 
                                                                                           "coords.x2"), c("min", "max")))
             , proj4string = new("CRS"
                                 , projargs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
             )
)

e <- extent(point@bbox)
e <- as(e,"SpatialPolygons")
e@proj4string <- point@proj4string

plot(e)
plot(point,add=T,col='red')

enter image description here

To extract what you want, use a loop selecting points by month. I created a list called plist to store each data.frame and also the month is added to each data.frame.

plist <- list()

for(i in unique(point$month)){
  p <- point[point$month == i,]
  index <- as.numeric(names(over(p,e)))
  plist[[i]] <- data.frame(as.data.frame(table(point$species[index])),month=i)
}

Check the resulting list with data.frames:

plist

## [[1]]
##   Var1 Freq month
## 1    1    2     1
## 2    2    1     1
## 3    3    1     1
## 
## [[2]]
##   Var1 Freq month
## 1    2    1     2
## 2    3    2     2
## 
## [[3]]
##   Var1 Freq month
## 1    1    1     3
## 2    2    1     3
## 3    3    1     3

Var1 is specie ID and Freq is count by specie.

0

Here is a solution based on dlply() from plyr embedded in a group_by-and-summarise pipeline from dplyr. Although I generally prefer to rely on fewer external packages in my code, I must admit that these two packages perform amazingly fast when it comes to such kinds of operations. I assume that your species data set has already been clipped, thus making an explicit call to over (or eg raster::crop) obsolete.

## sample data
library(sp)

dat = data.frame(long, lat, month, species)
coordinates(dat) = ~ long + lat # optionally set spatial context for future use
proj4string(dat) = "+init=epsg:4326"

## group into data.frames with species counts per month
library(plyr)
library(dplyr)

dat@data %>%
  group_by(month, species) %>%
  summarise(count = n()) %>%
  dlply("month")

# $`1`
# month species count
#     1       d     1
# 
# $`2`
# month species count
#     2       c     2
#     2       d     1
# 
# $`3`
# month species count
#     3       c     2
#     3       d     2
# 
# $`4`
# month species count
#     4       c     1

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.