How to get information, if plant species occur at specific coordinates by using shapefiles in R

Question

I want to create a tables using R that contain the information, if plant species occur at specific points (coordinates). Therefore, I want to use shapefiles, containing the distribution of several specific plant species. My output has to be a table, that for every point/coordinate indicate the presence of every plant species as 1 and the absence as 0.

First, I read in my first shapefile and the CVS table, containing the coordinates I need:

plant <- shapefile ('plant.shp')
birds<-read.csv2("bird_Coordinates.csv") 

Next, I extract the coordinates, save them in a dataframe and project the points over the distribution shapefile:

lats <- birds$lat
lons <- birds$lon
pts <- data.frame(x=lons,y=lats)
coordinates(pts) <- ~x+y
proj4string(pts) <- proj4string(plant)

When I now plot the shapefile and the coordinates, I see the shape of the plant distribution and two red points, which indicate that two of my about 60 points are within this distribution:

plot(plant)
points(pts, pch=20, col='red')

Next, I try to use over, to associate the points to the distribution. Here I used two different ways:

1.

over(pts, plant)$admin
cbind.data.frame(pts, plant=over(pts, plant)$admin)

Resulting in the warning message: Error in data.frame(..., check.names = FALSE) : Arguments imply differing numbers of rows: 64, 0

2.

plantsp <- !is.na(over(pts, as(plant, "SpatialPolygons")))
pts$plant <- over(pts, plant)$Unit_Name

Resulting the the warning message: Error in validObject(.Object) : invalid class “SpatialPointsDataFrame”: number of rows in data.frame and SpatialPoints don't match

So, both possibilities fail and I do not know what I did wrong. I know, that for this distribution range, only two points are inside the range, is the the source of the trouble? And how may I fix this? I would be so thankful, if anyone could tell me, how to get this cvs Table, containing the presence/absense information for the distribution ranges for every point!

Answer 1

I think this problem will be easier to solve with rgdal or sf (I recommend learning sf, as it's newer, but for now I'm more familiar with rgdal which is why I've used it here).

Load rgdal (which also loads sp, which we'll need):

library("rgdal")

# Open the shapefile and copy the projection,
# which we'll need for the bird data

plant <- readOGR(".", "N_columbiana")
wgs84 <- proj4string(plant)

# open the bird data and make a copy for when we change to
# spatialPointsDataFrame

birds <- read.csv2("bird_Coordinates.csv") 
birds_data <- birds

# Correct the order of long/lat so the coordinates are correct
birds <- birds[, c("lon", "lat")]
birds <- sp::SpatialPointsDataFrame(birds, birds_data,
                                    proj4string = CRS(wgs84))

# plot to make sure it's worked
plot(plant)
plot(birds, add = TRUE)

# the @data slot is now a data frame for the shapefile
# make all in_plant 0 by default before we subset
birds@data$in_plant <- 0

# Get the specimenID of those birds within the plant boundary
# the birds[plant, ] does the subsetting
# the @data$SpecimenID returns just the specimenID
in_plant <- birds[plant, ]@data$SpecimenID
birds@data$in_plant[in_plant] <- 1

# check it's all worked
plot(plant)
plot(birds[birds@data$in_plant == 1, ], add = TRUE)

If you want to extract just the data you can assign it to another object:

new_data <- birds@data
new_data
#     SpecimenID      lat         lon in_plant
#  1           1 46.06667    7.600000        0
#  2           2 46.60134    9.965973        0
#  3           3 46.02360    7.748607        0
#  4           4 46.60134    9.965973        0
#  5           5 46.91833   13.873611        0
#  ...
#  48         48 40.01861 -105.278056        1
#  49         49 40.02977 -105.581228        1
#  50         50 47.05917   13.615556        0