calculating area of most suitable raster habitat in R

Question

I have run Maxent for multiple species under present conditions and also under future climate change scenarios. I was quantifying changes between present and future suitable habitat using the nicheOverlap function and Schoener's D statistic. Quite a few of the organisms in my study are just moving farther up mountains so there is a lot of overlap as the future distribution is inside the present distribution (just occupying less area at higher elevations). By looking at the ascii files in QGIS I can see that there is less suitable habitat in terms of area in the future, so I want to quantify this. I have scoured the internet for a good way to calculate area for rasters and never found anything that perfectly suited my fancy. I therefore wrote up something that is an amalgamation of bits and pieces of various scripts. It is pasted below.

Two questions: 1) do you all agree this is doing what I think it is doing (calculating area in square kilometers)? 2) is there a way to simplify this? Specifically you'll see I go from a raster to a dataframe back to raster? Maybe I could stay in rasters?

Thanks for any input!

Rebecca

####
library(raster)

#load rasters
m <- raster("SpeciesA_avg.asc")
mf <- raster("SpeciesA_future_layers_avg.asc")

#change to dataframe
m.df <- as.data.frame(m, xy=TRUE)

#get rid of NAs
m.df1 <- na.omit(m.df)

#keep only cells that that have a suitability score above 0.5 (scores range from 0 to 1)
m.df2 <- m.df1[m.df1$SpeciesA_avg> 0.5,]

#re-rasterize just the suitable area
m.raster <- rasterFromXYZ(m.df2)

##same as above but for future projection
mf.df <- as.data.frame(mf, xy=TRUE)
mf.df1 <- na.omit(mf.df)
mf.df2 <- mf.df1[mf.df1$SpeciesA_future_layers_avg>0.5,]
mf.raster <-rasterFromXYZ(mf.df2)

#get sizes of all cells in current distribution raster
#note my original layers were 30 seconds or 1 km2. 
cell_size<-area(m.raster, na.rm=TRUE, weights=FALSE)

#delete NAs from all raster cells. It looks like these come back when switching from dataframe to raster
cell_size1<-cell_size[!is.na(cell_size)]

#compute area [km2] of all cells in raster
raster_area_present<-length(cell_size1)*median(cell_size1)
raster_area_present

#get sizes of all cells in future raster [km2]
cell_size<-area(mf.raster, na.rm=TRUE, weights=FALSE)

#delete NAs from vector of all raster cells
cell_size1<-cell_size[!is.na(cell_size)]

#compute area [km2] of all cells in geo_raster
raster_area_future<-length(cell_size1)*median(cell_size1)
raster_area_future

##calculate change in area
dif_area <- raster_area_present - raster_area_future
dif_area

Robert Hijmans Robert Hijmans · Accepted Answer · 2017-11-30T04:25:13

When you ask a question, you should provide a simple self-contained example. Not just dump your script that points to files we do not have. Writing a simple example teaches your R, and often helps you solve the problem by yourself. Anyway, I here is some example data and solution to your problem, I think:

library(raster)
#example data
m <- mf <- raster(ncol=10, nrow=10, vals=0)
m[,1] <- NA   
m[,3:7] <- 1
mf[,6:9] <- 1


# get rid of NAs (the example has none); should not be needed
m <- reclassify(m, cbind(NA, NA, 0))
mf <- reclassify(mf, cbind(NA, NA, 0))

# keep cells > 0.5 (scores range from 0 to 1)
m <- round(m)
mf <- round(mf)

# now combine the two layers, for example:
x <- m + mf * 10
# area of each cell
a <- area(x)
# sum area by class
z <- zonal(a, x, sum)

#     zone     value
#[1,]    0 152327547
#[2,]    1 152327547
#[3,]   10 101551698
#[4,]   11 101551698

zone 0 is "not current, nor future", 1 is "current only", 10 is "future only" and 11 is "current and future" The areas are in m^2.

You may want to check out this tutorial on maxent and other spatial distribution models: http://rspatial.org/sdm/

calculating area of most suitable raster habitat in R

1 Answers