发生建模后从栅格绘制多边形
Draw polygon from raster after occurrence modeling
我想使用 BIEN 使用的相同方法为物种出现绘制多边形，因此我可以同时使用我的多边形和他们的多边形。他们使用 Maxent 来模拟物种出现时的出现点。
因此，例如，这是一个 BIEN 多边形：
library(BIEN)
Mormolyca_ringens<- BIEN_ranges_load_species(species = "Mormolyca ringens")
#And this is a polygon, yes. A SpatialPolygonsDataFrame.
plot(wrld_simpl, xlim=c(-100,-40), ylim=c(-30,30), axes=TRUE,col="light yellow", bg="light blue")
plot(Mormolyca_ringens, col="green", add=TRUE)
Mormolyca ringens 多边形
好的，那我试着画我的多边形，因为 BIEN 缺少一些我需要的物种。
# first, you need to download the Maxent software here: http://biodiversityinformatics.amnh.org/open_source/maxent/
#and paste the "maxent.jar" file in the ’java’ folder of the ’dismo’ package, which is here:
system.file("java", package="dismo")
#You have to do this **before** loading the libraries
#install.packages("rJava")
library(rJava)
#If you get the message that cannot load this library, it's possible that your version of java is not 64bit.
#Go to Oracle and install Java for windows 64bit.
#If library still doesn't load: Look in your computer for the path where the java's jre file is and paste in the code below
Sys.setenv(JAVA_HOME="your\path\for\jre") #mine is "C:\Program Files\Java\jre1.8.0_144", for example
library(rJava)
library(dismo)
library(maptools)
#Giving credits: I wrote the following code based on this tutorial: https://cran.r-project.org/web/packages/dismo/vignettes/sdm.pdf
#Preparing the example data - the map
data(wrld_simpl)
ext = extent(-90, -32, -33, 23)
#Preparing the example data - presence data for Bradypus variegatus
file <- paste(system.file(package="dismo"), "/ex/bradypus.csv", sep="")
bradypus <- read.table(file, header=TRUE, sep=',')
bradypus <- bradypus[,-1] #don't need th first col
#Getting the predictors (the variables)
files <- list.files(path=paste(system.file(package="dismo"),
'/ex', sep=''), pattern='grd', full.names=TRUE )
predictors <- stack(files)
#making a training and a testing set.
group <- kfold(bradypus, 5)
pres_train <- bradypus[group != 1, ]
pres_test <- bradypus[group == 1, ]
#Creating the background
backg <- randomPoints(predictors, n=1000, ext=ext, extf = 1.25)
colnames(backg) = c('lon', 'lat')
group <- kfold(backg, 5)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]
# Running maxent
xm <- maxent(predictors, pres_train, factors='biome')
plot(xm)
#A response plot:
response(xm)
# Evaluating and predicting
e <- evaluate(pres_test, backg_test, xm, predictors)
px <- predict(predictors, xm, ext=ext, progress='text', overwrite=TRUE)
#Checking result of the prediction
par(mfrow=c(1,2))
plot(px, main='Maxent, raw values')
plot(wrld_simpl, add=TRUE, border='dark grey')
tr <- threshold(e, 'spec_sens')
plot(px > tr, main='presence/absence')
plot(wrld_simpl, add=TRUE, border='dark grey')
points(pres_train, pch='+')
此时，我有如下图像：
实例发生的预测
我正在尝试使用此代码从该栅格制作多边形：
predic_pol<-rasterToPolygons(px )
还有：
px_rec<-reclassify(px, rcl=0.5, include.lowest=FALSE)
px_pol<-rasterToPolygons(px_rec)
但我一直在获取我的范围的像素版本
能否请您给我一个提示，以便我可以从此光栅中提取多边形，就像 BIEN 的一样？ （我也是建模和 R 的新手......欢迎任何提示）
编辑：这是 px 控制台输出：
> px
class : RasterLayer
dimensions : 172, 176, 30272 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -120, -32, -56, 30 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : C:\Users\thai\Documents\ORCHIDACEAE\Ecologicos\w2\predictions\Trigonidiumobtusum_prediction.grd
names : layer
values : 6.705387e-06, 0.9999983 (min, max)
提前致谢
编辑 2：解决方案
感谢@Val，我做到了：
#Getting only the values>tr to make the polygon
#"tr" is what gives me the green raster instear of the multicolour one
pol <- rasterToPolygons(px>tr,function(x) x == 1,dissolve=T)
#Ploting
plot(wrld_simpl, xlim=c(-120,-20), ylim=c(-60,10), axes=TRUE,col="light yellow", bg="light blue")
plot(pol, add=T, col="green")
现在我得到了我想要的！谢谢！
（图中的多边形不一样只是因为我在得到@Val 的回答时使用了不同的数据集）
奖金问题：
你知道如何平滑边缘以便我得到一个非像素化的多边形吗？
我不了解 BIEN，所以我没有真正查看您示例的这一部分。我只是将您的 problem/question 概括为以下内容：
您有一个二进制栅格（0 表示不存在，1 表示存在）并且您想要将所有 1 的区域转换为多边形。
至于你的 px 光栅，你的值不是 0 和 1 而是更多 basically 0 和 basically 1. 但如果这是个问题，那可以很容易解决。
所以我试着用巴西地区重新创建你的例子：
library(raster)
library(rgeos)
# get Brasil borders
shp <- getData(country = 'BRA',level=0)
#create binary raster
r <- raster(extent(shp),resolution=c(0.5,0.5))
r[] <- NA # values have to be NA for the buffering
# take centroid of Brasil as center of species presence
cent <- gCentroid(shp)
# set to 1
r[cellFromXY(r,cent)] <- 1
# buffer presence
r <- buffer(r,width=1000000)
# set rest 0
r[is.na(r)] <- 0
# mask by borders
r <- mask(r,shp)
我猜这与您的光栅足够接近了：
所以现在要转换为多边形：
pol <- rasterToPolygons(r,function(x) x == 1,dissolve=T)
我使用一个函数来只获取值为 1 的像素。我还溶解了多边形，使其不是单个像素多边形，而是一个区域。有关其他选项，请参阅 rasterToPolygons。
现在一起绘制边界和新多边形：
plot(shp)
plot(pol,col='red',add=T)
这就是分布的多边形。这是控制台输出：
> pol
class : SpatialPolygonsDataFrame
features : 1
extent : -62.98971, -43.48971, -20.23512, -1.735122 (xmin, xmax, ymin, ymax)
coord. ref. : NA
variables : 1
names : layer
min values : 1
max values : 1
希望对您有所帮助！
编辑：奖励答案
您必须清楚，多边形的像素化边界代表了数据的准确表示。所以任何改变都意味着精度的损失。现在，根据您的目的，这可能无关紧要。
有多种方法可以实现它，可以是在栅格端进行分解和 smoothing/filtering 等，也可以是在多边形端，您可以在其中对多边形应用特定的过滤器，例如 this。
如果纯粹是为了美观，你可以尝试 rgeos 包中的 gSimplify：
# adjust tol for smoothness
pol_sm <- gSimplify(pol,tol=0.5)
plot(pol)
lines(pol_sm,col='red',lwd=2)
我想使用 BIEN 使用的相同方法为物种出现绘制多边形，因此我可以同时使用我的多边形和他们的多边形。他们使用 Maxent 来模拟物种出现时的出现点。
因此，例如，这是一个 BIEN 多边形：
library(BIEN)
Mormolyca_ringens<- BIEN_ranges_load_species(species = "Mormolyca ringens")
#And this is a polygon, yes. A SpatialPolygonsDataFrame.
plot(wrld_simpl, xlim=c(-100,-40), ylim=c(-30,30), axes=TRUE,col="light yellow", bg="light blue")
plot(Mormolyca_ringens, col="green", add=TRUE)
Mormolyca ringens 多边形
好的，那我试着画我的多边形，因为 BIEN 缺少一些我需要的物种。
# first, you need to download the Maxent software here: http://biodiversityinformatics.amnh.org/open_source/maxent/
#and paste the "maxent.jar" file in the ’java’ folder of the ’dismo’ package, which is here:
system.file("java", package="dismo")
#You have to do this **before** loading the libraries
#install.packages("rJava")
library(rJava)
#If you get the message that cannot load this library, it's possible that your version of java is not 64bit.
#Go to Oracle and install Java for windows 64bit.
#If library still doesn't load: Look in your computer for the path where the java's jre file is and paste in the code below
Sys.setenv(JAVA_HOME="your\path\for\jre") #mine is "C:\Program Files\Java\jre1.8.0_144", for example
library(rJava)
library(dismo)
library(maptools)
#Giving credits: I wrote the following code based on this tutorial: https://cran.r-project.org/web/packages/dismo/vignettes/sdm.pdf
#Preparing the example data - the map
data(wrld_simpl)
ext = extent(-90, -32, -33, 23)
#Preparing the example data - presence data for Bradypus variegatus
file <- paste(system.file(package="dismo"), "/ex/bradypus.csv", sep="")
bradypus <- read.table(file, header=TRUE, sep=',')
bradypus <- bradypus[,-1] #don't need th first col
#Getting the predictors (the variables)
files <- list.files(path=paste(system.file(package="dismo"),
'/ex', sep=''), pattern='grd', full.names=TRUE )
predictors <- stack(files)
#making a training and a testing set.
group <- kfold(bradypus, 5)
pres_train <- bradypus[group != 1, ]
pres_test <- bradypus[group == 1, ]
#Creating the background
backg <- randomPoints(predictors, n=1000, ext=ext, extf = 1.25)
colnames(backg) = c('lon', 'lat')
group <- kfold(backg, 5)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]
# Running maxent
xm <- maxent(predictors, pres_train, factors='biome')
plot(xm)
#A response plot:
response(xm)
# Evaluating and predicting
e <- evaluate(pres_test, backg_test, xm, predictors)
px <- predict(predictors, xm, ext=ext, progress='text', overwrite=TRUE)
#Checking result of the prediction
par(mfrow=c(1,2))
plot(px, main='Maxent, raw values')
plot(wrld_simpl, add=TRUE, border='dark grey')
tr <- threshold(e, 'spec_sens')
plot(px > tr, main='presence/absence')
plot(wrld_simpl, add=TRUE, border='dark grey')
points(pres_train, pch='+')
此时，我有如下图像：
实例发生的预测
我正在尝试使用此代码从该栅格制作多边形：
predic_pol<-rasterToPolygons(px )
还有：
px_rec<-reclassify(px, rcl=0.5, include.lowest=FALSE)
px_pol<-rasterToPolygons(px_rec)
但我一直在获取我的范围的像素版本
能否请您给我一个提示，以便我可以从此光栅中提取多边形，就像 BIEN 的一样？ （我也是建模和 R 的新手......欢迎任何提示）
编辑：这是 px 控制台输出：
> px
class : RasterLayer
dimensions : 172, 176, 30272 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -120, -32, -56, 30 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : C:\Users\thai\Documents\ORCHIDACEAE\Ecologicos\w2\predictions\Trigonidiumobtusum_prediction.grd
names : layer
values : 6.705387e-06, 0.9999983 (min, max)
提前致谢
编辑 2：解决方案
感谢@Val，我做到了：
#Getting only the values>tr to make the polygon
#"tr" is what gives me the green raster instear of the multicolour one
pol <- rasterToPolygons(px>tr,function(x) x == 1,dissolve=T)
#Ploting
plot(wrld_simpl, xlim=c(-120,-20), ylim=c(-60,10), axes=TRUE,col="light yellow", bg="light blue")
plot(pol, add=T, col="green")
现在我得到了我想要的！谢谢！ （图中的多边形不一样只是因为我在得到@Val 的回答时使用了不同的数据集）
奖金问题：
你知道如何平滑边缘以便我得到一个非像素化的多边形吗？
我不了解 BIEN，所以我没有真正查看您示例的这一部分。我只是将您的 problem/question 概括为以下内容：
您有一个二进制栅格（0 表示不存在，1 表示存在）并且您想要将所有 1 的区域转换为多边形。
至于你的 px 光栅，你的值不是 0 和 1 而是更多 basically 0 和 basically 1. 但如果这是个问题，那可以很容易解决。
所以我试着用巴西地区重新创建你的例子：
library(raster)
library(rgeos)
# get Brasil borders
shp <- getData(country = 'BRA',level=0)
#create binary raster
r <- raster(extent(shp),resolution=c(0.5,0.5))
r[] <- NA # values have to be NA for the buffering
# take centroid of Brasil as center of species presence
cent <- gCentroid(shp)
# set to 1
r[cellFromXY(r,cent)] <- 1
# buffer presence
r <- buffer(r,width=1000000)
# set rest 0
r[is.na(r)] <- 0
# mask by borders
r <- mask(r,shp)
我猜这与您的光栅足够接近了：
所以现在要转换为多边形：
pol <- rasterToPolygons(r,function(x) x == 1,dissolve=T)
我使用一个函数来只获取值为 1 的像素。我还溶解了多边形，使其不是单个像素多边形，而是一个区域。有关其他选项，请参阅 rasterToPolygons。
现在一起绘制边界和新多边形：
plot(shp)
plot(pol,col='red',add=T)
这就是分布的多边形。这是控制台输出：
> pol
class : SpatialPolygonsDataFrame
features : 1
extent : -62.98971, -43.48971, -20.23512, -1.735122 (xmin, xmax, ymin, ymax)
coord. ref. : NA
variables : 1
names : layer
min values : 1
max values : 1
希望对您有所帮助！
编辑：奖励答案
您必须清楚，多边形的像素化边界代表了数据的准确表示。所以任何改变都意味着精度的损失。现在，根据您的目的，这可能无关紧要。
有多种方法可以实现它，可以是在栅格端进行分解和 smoothing/filtering 等，也可以是在多边形端，您可以在其中对多边形应用特定的过滤器，例如 this。
如果纯粹是为了美观，你可以尝试 rgeos 包中的 gSimplify：
# adjust tol for smoothness
pol_sm <- gSimplify(pol,tol=0.5)
plot(pol)
lines(pol_sm,col='red',lwd=2)