diff --git a/DESCRIPTION b/DESCRIPTION
index 2feeb5f..bc4d1d7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -45,7 +45,9 @@ Imports:
     spatstat.random,
     spatstat.geom,
     stats,
-    tibble
+    tibble,
+    Rcpp,
+    terra
 Suggests:
     ggplot2,
     highcharter,
diff --git a/R/nlm_percolation.R b/R/nlm_percolation.R
index 25f27b5..294f551 100755
--- a/R/nlm_percolation.R
+++ b/R/nlm_percolation.R
@@ -10,6 +10,8 @@
 #' Resolution of the raster.
 #' @param prob [\code{numerical(1)}]\cr
 #' Probability value for setting a cell to 1.
+#' @param crs [\code{character(1)}]\cr
+#' Coordinate reference system for new raster. If blank, defaults to WGS 84.
 #'
 #' @details
 #' The simulation of a random percolation map is accomplished in two steps:
@@ -22,7 +24,7 @@
 #'  TRUE - if it is higher the cell is set to FALSE.}
 #' }
 #'
-#' @return RasterLayer
+#' @return SpatRaster
 #'
 #' @examples
 #' # simulate percolation model
@@ -49,7 +51,8 @@
 nlm_percolation <- function(ncol,
                             nrow,
                             resolution = 1,
-                            prob = 0.5) {
+                            prob = 0.5,
+                            crs = "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") {
 
   # Check function arguments ----
   checkmate::assert_count(ncol, positive = TRUE)
@@ -57,27 +60,16 @@ nlm_percolation <- function(ncol,
   checkmate::assert_numeric(resolution)
   checkmate::assert_true(prob <= 1, na.ok = FALSE)
   checkmate::assert_true(prob >= 0, na.ok = FALSE)
-
-  percolation_matrix <- matrix(NA, nrow = nrow, ncol = ncol)
-
-  percolation_matrix[] <- vapply(
-    percolation_matrix,
-    function(x) {
-      ifelse(stats::runif(1, 0, 1) < prob, TRUE, FALSE)
-    },
-    logical(1)
-  )
-
-  percolation_raster <-
-    raster::raster(percolation_matrix)
-
-  # specify resolution ----
-  raster::extent(percolation_raster) <- c(
-    0,
-    ncol(percolation_raster) * resolution,
-    0,
-    nrow(percolation_raster) * resolution
-  )
+    
+  percolation_raster = terra::rast(nrows = nrow, ncols = ncol, nlyrs=1,
+                                   resolution = resolution,
+                                   extent = c(0, ncol * resolution,
+                                              0, nrow * resolution),
+                                   crs = crs,
+                                   vals = sample(c(0,1), 
+                                                 size = ncol * nrow, 
+                                                 replace = T, 
+                                                 prob = c(1-prob, prob)))
 
   return(percolation_raster)
 }
diff --git a/R/nlm_randomcluster.R b/R/nlm_randomcluster.R
index a1ce515..73d7f9e 100644
--- a/R/nlm_randomcluster.R
+++ b/R/nlm_randomcluster.R
@@ -18,8 +18,10 @@
 #' This directly controls the number of types via the given length.
 #' @param rescale [\code{logical(1)}]\cr
 #' If \code{TRUE} (default), the values are rescaled between 0-1.
+#' @param crs [\code{character(1)}]\cr
+#' The crs for the new raster. If blank, defaults to WGS 84.
 #'
-#' @return Raster with random values ranging from 0-1.
+#' @return SpatRaster with random values ranging from 0-1.
 #'
 #' @details
 #' This is a direct implementation of steps A - D of the modified random clusters algorithm
@@ -52,7 +54,8 @@ nlm_randomcluster <- function(ncol, nrow,
                               p,
                               ai = c(0.5, 0.5),
                               neighbourhood = 4,
-                              rescale = TRUE) {
+                              rescale = TRUE, 
+                              crs = "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") {
 
   # Check function arguments ----
   checkmate::assert_count(ncol, positive = TRUE)
@@ -65,26 +68,29 @@ nlm_randomcluster <- function(ncol, nrow,
   checkmate::assert_logical(rescale)
 
   # Step A - Create percolation map
-  ranclumap <- nlm_percolation(ncol, nrow, p, resolution = resolution)
+  ranclumap <- nlm_percolation(ncol, nrow, p, resolution = resolution, crs = crs)
 
   # Step B - Cluster identification (clustering of adjoining pixels)
-  ranclumap <- raster::clump(ranclumap, direction = neighbourhood, gaps = FALSE)
-
+  ranclumap <- terra::patches(x = ranclumap, directions = neighbourhood, zeroAsNA = T)
+  #terra::values(ranclumap)[is.na(terra::values(ranclumap))] <- 0
+  
   # Step C - Cluster type assignation
   # number of different cluster
-  numclu <- max(raster::values(ranclumap), na.rm = TRUE)
+  numclu <- length(unique(ranclumap[][!is.na(ranclumap[])]))
   # assign to each cluster nr a new category given by Ai
   clutyp <- sample(seq_along(ai), numclu, replace = TRUE, prob = ai)
   # write back new category nr
-  raster::values(ranclumap) <- clutyp[raster::values(ranclumap)]
+  class_df <- data.frame('is' = unique(ranclumap[][!is.na(ranclumap[])]),
+                         'becomes' = clutyp)
+  ranclumap <- terra::classify(ranclumap, rcl = class_df, include.lowest=T)
 
   # Step D - Filling the map
   # helperfuction to choose values
   fillit <- function(cid) {
     # get neighbour cells
-    nbrs <- raster::adjacent(ranclumap, cid, directions = 8, pairs = FALSE)
+    nbrs <- terra::adjacent(ranclumap, cid, directions = 8, pairs = FALSE)
     # count neighbour values (exclude NA see Saura 2000 paper)
-    vals <- table(raster::values(ranclumap)[nbrs])
+    vals <- table(terra::values(ranclumap)[nbrs])
     # check if everything is NA
     if (!length(vals)) {
       # be a rebel get your own value
@@ -103,24 +109,24 @@ nlm_randomcluster <- function(ncol, nrow,
 
   # identify unfilled cells
   gaps <- dplyr::rowwise(tibble::tibble(
-    ctf = (1:(ncol * nrow))[is.na(raster::values(ranclumap))]
+    ctf = (1:(ncol * nrow))[is.na(terra::values(ranclumap))]
     ))
   # get values for the gaps
   gaps <- dplyr::mutate(gaps, val = fillit(ctf))
   # feed it back in the map
-  raster::values(ranclumap)[gaps$ctf] <- gaps$val
+  terra::values(ranclumap)[gaps$ctf] <- gaps$val
 
   # specify resolution ----
-  raster::extent(ranclumap) <- c(
-    0,
-    ncol(ranclumap) * resolution,
-    0,
-    nrow(ranclumap) * resolution
-  )
+  # raster::extent(ranclumap) <- c(
+  #   0,
+  #   ncol(ranclumap) * resolution,
+  #   0,
+  #   nrow(ranclumap) * resolution
+  # )
 
   # Rescale values to 0-1 ----
   if (rescale == TRUE) {
-    ranclumap <- util_rescale(ranclumap)
+    ranclumap <- ranclumap / max(ranclumap[])
   }
 
   return(ranclumap)
diff --git a/tests/testthat.R b/tests/testthat.R
index 3c66efe..d561581 100644
--- a/tests/testthat.R
+++ b/tests/testthat.R
@@ -1,4 +1,4 @@
-library(testthat)
-library(NLMR)
-
-test_check("NLMR")
+library(testthat)
+library(NLMR)
+
+test_check("NLMR")
diff --git a/tests/testthat/test_percolation.R b/tests/testthat/test_percolation.R
index 0c245bd..94901af 100644
--- a/tests/testthat/test_percolation.R
+++ b/tests/testthat/test_percolation.R
@@ -1,22 +1,22 @@
-# nolint start
-context("nlm_percolation")
-
-percolation  <- nlm_percolation(ncol = 9, nrow = 12, prob = 0.5)
-
-test_that("nlm_percolation behaves like it should", {
-  expect_that(percolation , is_a("RasterLayer"))
-})
-
-test_that("nlm_percolation produces the right number of rows", {
-  expect_equal(percolation@nrows, 12)
-})
-
-test_that("nlm_percolation produces the right number of columns", {
-  expect_equal(percolation@ncols, 9)
-})
-
-test_that("nlm_percolation produces the right number of columns", {
-  expect_equal(length(unique(percolation@data@values)), 2)
-})
-
-# nolint end
+# nolint start
+context("nlm_percolation")
+
+percolation  <- nlm_percolation(ncol = 9, nrow = 12, prob = 0.5)
+
+test_that("nlm_percolation behaves like it should", {
+  expect_that(percolation , is_a("SpatRaster"))
+})
+
+test_that("nlm_percolation produces the right number of rows", {
+  expect_equal(terra::nrow(percolation), 12)
+})
+
+test_that("nlm_percolation produces the right number of columns", {
+  expect_equal(terra::ncol(percolation), 9)
+})
+
+test_that("nlm_percolation produces the right number of columns", {
+  expect_equal(length(unique(terra::values(percolation))), 2)
+})
+
+# nolint end
diff --git a/tests/testthat/test_randomcluster.R b/tests/testthat/test_randomcluster.R
index cc2dcba..5f932bd 100644
--- a/tests/testthat/test_randomcluster.R
+++ b/tests/testthat/test_randomcluster.R
@@ -1,19 +1,38 @@
-# nolint start
-context("nlm_randomcluster")
-
-suppressWarnings(random_cluster  <- nlm_randomcluster(ncol = 40, nrow = 30,
-                                     neighbourhood = 4, p = 0.4))
-
-test_that("nlm_randomcluster behaves like it should", {
-  expect_that(random_cluster , is_a("RasterLayer"))
-})
-
-test_that("nlm_randomcluster produces the right number of rows", {
-  expect_equal(random_cluster@nrows, 30)
-})
-
-test_that("nlm_randomcluster produces the right number of columns", {
-  expect_equal(random_cluster@ncols, 40)
-})
-
-# nolint end
+# nolint start
+context("nlm_randomcluster")
+
+suppressWarnings(random_cluster  <- nlm_randomcluster(ncol = 300, nrow = 250,
+                                     neighbourhood = 4, p = 0.4, ai = c(0.1, 0.3, 0.6),
+                                     rescale = F))
+
+test_that("nlm_randomcluster behaves like it should", {
+  expect_that(random_cluster , is_a("SpatRaster"))
+})
+
+test_that("nlm_randomcluster produces the right number of rows", {
+  expect_equal(terra::nrow(random_cluster), 250)
+})
+
+test_that("nlm_randomcluster produces the right number of columns", {
+  expect_equal(terra::ncol(random_cluster), 300)
+})
+
+test_that("nlm_randomcluster produces expected values", {
+  expect_equal(length(unique(terra::values(random_cluster))), 3)
+})
+
+test_that("nlm_randomcluster produces proportions within 0.05 of expected", {
+  expect_equal(terra::freq(random_cluster)[,3][1] / length(terra::values(random_cluster)), 
+               0.1, tolerance = 0.05)
+})
+
+test_that("nlm_randomcluster produces proportions within 0.05 of expected", {
+  expect_equal(terra::freq(random_cluster)[,3][2] / length(terra::values(random_cluster)), 
+               0.3, tolerance = 0.05)
+})
+
+test_that("nlm_randomcluster produces proportions within 0.05 of expected", {
+  expect_equal(terra::freq(random_cluster)[,3][3] / length(terra::values(random_cluster)), 
+               0.6, tolerance = 0.05)
+})
+# nolint end