计算R中最合适的栅格栖息地的面积
calculating area of most suitable raster habitat in R
我有 运行 当前条件下以及未来气候变化情景下多个物种的 Maxent。我正在使用 nicheOverlap 函数和 Schoener 的 D 统计量来量化当前和未来合适栖息地之间的变化。在我的研究中，有相当多的生物体只是向山上移动得更远，所以有很多重叠，因为未来的分布在现在的分布范围内（只是在更高的海拔占据更少的面积）。通过查看 QGIS 中的 ascii 文件，我可以看到未来在面积方面的适宜栖息地较少，因此我想对其进行量化。我在互联网上搜索了一种计算栅格面积的好方法，但从未找到任何完全适合我的想法。因此，我写了一些东西，它是各种脚本的点点滴滴的融合。下面贴出来了。
两个问题：
1）你们都同意这是在做我认为正在做的事情（以平方公里计算面积）吗？
2）有没有办法简化这个？具体来说，您会看到我从栅格到数据框再回到栅格？也许我可以留在光栅中？
感谢任何意见！
丽贝卡
####
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
当你提出问题时，你应该提供一个简单的自包含示例。不只是转储指向我们没有的文件的脚本。编写一个简单的示例可以教会您 R，并且通常可以帮助您自己解决问题。无论如何，我认为这是一些示例数据和解决您的问题的方法：
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
区域 0 是 "not current, nor future"，区域 1 是 "current only"，区域 10 是 "future only"，区域 11 是 "current and future"
面积以 m^2 为单位。
您可能想查看有关 maxent 和其他空间分布模型的教程：http://rspatial.org/sdm/
我有 运行 当前条件下以及未来气候变化情景下多个物种的 Maxent。我正在使用 nicheOverlap 函数和 Schoener 的 D 统计量来量化当前和未来合适栖息地之间的变化。在我的研究中，有相当多的生物体只是向山上移动得更远，所以有很多重叠，因为未来的分布在现在的分布范围内（只是在更高的海拔占据更少的面积）。通过查看 QGIS 中的 ascii 文件，我可以看到未来在面积方面的适宜栖息地较少，因此我想对其进行量化。我在互联网上搜索了一种计算栅格面积的好方法，但从未找到任何完全适合我的想法。因此，我写了一些东西，它是各种脚本的点点滴滴的融合。下面贴出来了。
两个问题： 1）你们都同意这是在做我认为正在做的事情（以平方公里计算面积）吗？ 2）有没有办法简化这个？具体来说，您会看到我从栅格到数据框再回到栅格？也许我可以留在光栅中？
感谢任何意见！
丽贝卡
####
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
当你提出问题时，你应该提供一个简单的自包含示例。不只是转储指向我们没有的文件的脚本。编写一个简单的示例可以教会您 R，并且通常可以帮助您自己解决问题。无论如何，我认为这是一些示例数据和解决您的问题的方法：
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
区域 0 是 "not current, nor future"，区域 1 是 "current only"，区域 10 是 "future only"，区域 11 是 "current and future" 面积以 m^2 为单位。
您可能想查看有关 maxent 和其他空间分布模型的教程：http://rspatial.org/sdm/