I want to draw polygons for species occurrence using the same methods BIEN uses, so I can use both my polygons and theirs. They use Maxent to model species occurrence when they have more then occurrence points.

So, this is, for example, a BIEN polygon:

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 polygon

Ok, then I'm trying to draw my polygons because BIEN lacks some for species I need.

# 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='+')

At this point, I have the following image:

Prediction for example's occurrence

And I'm trying to make a polygon from this raster with this code:

predic_pol<-rasterToPolygons(px )

And also:

px_rec<-reclassify(px, rcl=0.5,  include.lowest=FALSE)
px_pol<-rasterToPolygons(px_rec)

But i keep getting a pixels version of my extent

Can you please give me a hint so I can extract a polygon out of this raster, like the BIEN's one? (Also I'm new to modeling and to R... any tips are welcome)

EDIT: this is the px console output:

> 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)

Thank you in advance

Edit 2: Solution

Thanks to @Val I got to this:

#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")

And now I have what I wanted! Thank you! (The polygon is not the same in the figures only because I used a different data set I had at my environment at the moment I got @Val 's answer)

Bonus question:

Do you know how to smooth the edges so I get a non pixelized polygon?

I don't know BIEN, so I din't really look at this part of your example. I just generalized your problem/question down to the following:

You have a binary raster (with 0 for absence and 1 for presence) and you want to convert all areas with 1 to a polygon.

As for your px raster, it's a bit odd that your values are not 0 and 1 but more basically 0 and basically 1. But if that's a problem, that can be an easy fix.

So I tried to recreate your example with just the area of Brasil:

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)

This is close enough to your raster I guess:

So now to the conversion to the polygon:

pol <- rasterToPolygons(r,function(x) x == 1,dissolve=T)

I use a function to only get pixels with value 1. Also I dissolve the polygons to not have single pixel polygons but rather an area. See rasterToPolygons for other options.

And now plot the borders and the new polygon together:

plot(shp)
plot(pol,col='red',add=T)

And there you have it, a polygon of the distribution. This is the console output:

> 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 

Hope that helps!

Edit: Bonus answer

You have to be clear, that the pixelized boundaries of your polygon(s) represent an accurate representation of your data. So any change to that means a loss of precision. Now, depending on your purpose, that might not matter.

There's multiple ways to achieve it, either at the raster side with disaggregating and smoothing/filtering etc. or at the polygon side, where you can apply specific filters to the polygons like this.

If it's purely aesthetic, you can try gSimplify from the rgeos package:

# adjust tol for smoothness
pol_sm <- gSimplify(pol,tol=0.5)

plot(pol)
lines(pol_sm,col='red',lwd=2)