What you need is to have only one xy for all points which occur in the same grid cell.
What I did was:
- create rasters with the distribution of each species
- load all rasters and stack them
- get minimum value from the stack (it will get only values from cells that have at least 1 sp present)
- transform this raster with min value to xy coordinates
- use xy coordinates to extract from raster
Let's take it step by step:
1. rasterize distribution of all species:
What you need:
- sp.xy - a spatial points data frame, with at least 4 columns: x, y, Taxon, Presence (corresponding to long, lat, species name, and value "1")
- inside sp.xy you need this column named Presence, with value 1 for all species. It is very easy to add if you start with an xy excel file (species name, long, lat - you just add another column with value 1 and name Presence).
- you can transform this excel file to an spdf file (spatial) using the latter part of this code (I did it for an xy file here; if you have problems with it, just search how to do that, you will find help)
- You also need a raster with the resolution and extent of the raster layers you wish to extract values from. You can actually load one of these rasters in R, and use it. It will not change it on the disk drive. It is referenced in my code with the name "raster". You can call it that or whaterver, just know what it needs to be.
Rasterize species data
get species names
sp.names=unique(sp.xy$Taxon)
rasterize each species using the presence field, and the function "min" in the rasterize function:
for(i in 1:length(sp.names))
{
tmp1=sp.xy[sp.xy$Taxon == sp.names[[i]],]
tmp2=rasterize(tmp1, raster, field="Presence", fun="min")
out.name=sp.names[[i]]
setwd("C:/output") ## insert your folder path here
writeRaster(tmp2, filename=as.character(out.name), format="GTiff", overwrite=TRUE)
}
2. Load all outputed rasters
require(raster)
setwd("C:/output") ## edit here with your path
names=list.files("C:/output", pattern="tif$", full.names=FALSE)
for (tif in names) assign(tif, raster(tif))
spp=stack(mget(names))
3. get xy of center of raster cells with value "1"... this can be done manually for each species or in a loop (got lazy here, you have an example for 1 species)
pts.sp1=rasterToPoints(names[[1]], fun=function(x){x==1})
4. transform to spdf (WGS - replace with proper CRS if you need)
pts.sp1=as.data.frame(pts.sp1)
coordinates(pts.sp1) <- ~x+y
projection(pts.sp1) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
5. extract from your raster stack
var.values=extract(raster, pts.sp1)