3

I have created several raster layers within ArcMap that I wish to extract values from in order to analyse within R.

However, I'm keen to avoid pseudo replication of my species presence data and therefore only wish to have one point value per cell.

Is there a way to do this?

2 Answers 2

5

What you need is to have only one xy for all points which occur in the same grid cell.

What I did was:

  1. create rasters with the distribution of each species
  2. load all rasters and stack them
  3. get minimum value from the stack (it will get only values from cells that have at least 1 sp present)
  4. transform this raster with min value to xy coordinates
  5. 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)
1

"Extract value to point" will extract the cell values for each point. If you want to remove duplicate points (several points in one cell), you could compute the row and column values of each point in the coordiante system of the grid cell, then you can identify duplicate entries (using "Find identical").

for the corresponding row/column values, take the X (Y) coordinates of your points, remove the left (bottom) coordinate of your grid, then divide by the cell size (integer division). If points have the same X and Y, they are in the same same. Watch out that some raster formats have their corners measured at the center of the pixel, so make sure that you check the coordinates.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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