diff --git a/NAMESPACE b/NAMESPACE
index df694c3..5dcfeaf 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -8,6 +8,7 @@ export(nm_eval)
 export(nm_gbif_suggestion)
 export(nm_find_hcv)
 export(nm_predict)
+export(nm_extract)
 ## Import
 importFrom("grDevices", "gray.colors", "rainbow")
 importFrom("graphics", "points")
diff --git a/R/nm_extract.R b/R/nm_extract.R
new file mode 100644
index 0000000..b7579ea
--- /dev/null
+++ b/R/nm_extract.R
@@ -0,0 +1,72 @@
+#' Extract raster data values based on point coordinates (longitude, latitude)
+#'
+#' @param data data.frame. Database with species presence, presence-absence,
+#' or pseudo-absence records with longitude and latitude coordinates.
+#' @param longitude character. Column name with spatial longitude coordinates.
+#' @param latitude character. Column name with spatial latitude coordinates.
+#' @param env_layer SpatRaster. Raster or raster stack.
+#' @param current_crs character. Current coordinate reference system of coordinates. One of:
+#' * `geographic`: Apllicable for longitude and latitude in decimal
+#' (e.g 31.796406, -24.536200)
+#' * `projected`: Applicable for longitude and latitude in projected
+#' system (e.g 7280635, 982148 - UTM zone 35S)
+#' @param variables character. Vector with the raster names to extract value from.
+#' If NULL (default) the function will return data for all raster.
+#' @param filter_na logical. If filter_na = TRUE (default), the rows with NA values
+#' for any of the raster are removed from the returned tibble.
+#'
+#' @return tibble with original data bound with extracted data
+#'
+#'
+#' @examples
+#'  library(terra)
+#'  wbv_path <- paste0(system.file("extdata", package = "nimo"), "/WbV_subset_occ_KNP.csv")
+#'  wbv_df <- read.csv(wbv_path)
+#'
+#'  env_layers_path <- paste0(system.file("extdata", package = "nimo"), "/env_layers")
+#'  env_layers <- terra::rast(list.files(env_layers_path, full.names = TRUE))
+#'
+#' extracted <- nm_extract(data = wbv_df, longitude = "decimalLongitude",
+#' latitude = "decimalLatitude", env_layer = env_layers)
+#'
+#' @export
+nm_extract <- function (data,
+                        longitude,
+                        latitude,
+                        env_layer,
+                        current_crs = "geographic",
+                        variables = NULL,
+                        filter_na = TRUE)
+{
+  if (is.null(variables)) {
+    variables <- names(env_layer)
+  }
+  data <- data %>%
+    dplyr::filter(!is.na(!!dplyr::sym(longitude))) %>%
+    dplyr::filter(!is.na(!!dplyr::sym(latitude)))
+
+  if (! current_crs %in% c("geographic", "projected")) {
+    stop("The current crs must be one of 'geographic' and 'projected'")
+  }
+  if(current_crs == "geographic"){
+    coord_rs <- 4326
+  }else{
+    coord_rs <- terra::crs(env_layer)
+  }
+  sp_data <- sf::st_as_sf(data, coords = c(longitude, latitude), crs = coord_rs) %>%
+    sf::st_transform(crs = terra::crs(env_layer))
+  extract_data <- dplyr::tibble(data,
+                                terra::extract(env_layer[[variables]], sp_data, cells = FALSE) %>%
+                                  dplyr::select({{variables}}))
+  if (filter_na) {
+    complete_vec <- stats::complete.cases(extract_data[,
+                                                       variables])
+    if (sum(!complete_vec) > 0) {
+      message(sum(!complete_vec), " rows were excluded from database because NAs were found")
+      extract_data <- extract_data %>% dplyr::filter(complete_vec)
+    }
+  }
+  return(extract_data)
+}
+
+
diff --git a/inst/nimo/server.R b/inst/nimo/server.R
index 3ee746c..f0e782c 100644
--- a/inst/nimo/server.R
+++ b/inst/nimo/server.R
@@ -437,39 +437,41 @@ occ_dt <- reactive({
 
   reduce_colin <- eventReactive(input$reduce_collin, {
     req(env_layers())
-    if(input$coli_method == "pearson"){
-      colin_var <- correct_colinvar(env_layer = env_layers(),
-                                    method = c("pearson", th = as.character(input$pearson_threshold)))
-      cr_df <- colin_var$cor_table # table
-      enlayer <- colin_var$cor_variables # layer
-      highly_corr_vars <- nm_find_hcv(cr_df, cutoff = input$pearson_threshold)
-      rm_enlayer <- highly_corr_vars
-    } else if(input$coli_method == "vif"){
-      colin_var <- correct_colinvar(env_layer = env_layers(), method = c("vif", th = as.character(input$vif_threshold)))
-      cr_df <- colin_var$vif_table # table
-      enlayer <- list()
-      for (i in 1:terra::nlyr(colin_var$env_layer)) {
-        l <- colin_var$env_layer[[i]]
-        enlayer[[i]] <- l
-      }
-      rm_enlayer <- colin_var$removed_variables
-    } else if(input$coli_method == "pca"){
-      colin_var <- correct_colinvar(env_layer = env_layers(), method = c("pca"))
-      cr_df <- colin_var$coefficients # table
-      enlayer <- list()
-      for (i in 1:terra::nlyr(colin_var$env_layer)) {
-        l <- colin_var$env_layer[[i]]
-        enlayer[[i]] <- l
-      }
-      rm_enlayer <- colin_var$cumulative_variance
-    } else if(input$coli_method == "fa" & terra::nlyr(env_layers()) > 2){
+    tryCatch({
+      if(input$coli_method == "pearson"){
+        colin_var <- correct_colinvar(env_layer = env_layers(),
+                                      method = c("pearson", th = as.character(input$pearson_threshold)))
+        cr_df <- colin_var$cor_table # table
+        enlayer <- colin_var$cor_variables # layer
+        highly_corr_vars <- nm_find_hcv(cr_df, cutoff = input$pearson_threshold)
+        rm_enlayer <- highly_corr_vars
+      } else if(input$coli_method == "vif"){
+        colin_var <- correct_colinvar(env_layer = env_layers(), method = c("vif", th = as.character(input$vif_threshold)))
+        cr_df <- colin_var$vif_table # table
+        enlayer <- list()
+        for (i in 1:terra::nlyr(colin_var$env_layer)) {
+          l <- colin_var$env_layer[[i]]
+          enlayer[[i]] <- l
+        }
+        rm_enlayer <- colin_var$removed_variables
+      } else if(input$coli_method == "pca"){
+        colin_var <- correct_colinvar(env_layer = env_layers(), method = c("pca"))
+        cr_df <- colin_var$coefficients # table
+        enlayer <- list()
+        for (i in 1:terra::nlyr(colin_var$env_layer)) {
+          l <- colin_var$env_layer[[i]]
+          enlayer[[i]] <- l
+        }
+        rm_enlayer <- colin_var$cumulative_variance
+      } else if(input$coli_method == "fa" & terra::nlyr(env_layers()) > 2){
         colin_var <- correct_colinvar(env_layer = env_layers(), method = c("fa"))
         cr_df <- colin_var$loadings # table
         enlayer <- colin_var$env_layer
         rm_enlayer <- colin_var$removed_variables
-    }
+      }
 
       return(list(cr_df, enlayer, rm_enlayer, colin_var))
+    }, error = error)
   })
 
   cr_env <- eventReactive(input$reduce_collin, {
@@ -955,33 +957,30 @@ observeEvent(input$valided_dp, {
   ## update predictor selection
 kept_colin_var <- reactiveValues(rc = NULL)
 
-rc <- reactive({
-  req(env_layers())
-  kept_colin_var$rc <- reduce_colin()[[3]]
-})
-
-
   observe({
     req(env_layers())
-    print(kept_colin_var$rc)
-    #print(kept_colin_var)
-    if (is.null(rc())) {
-      reduce_colin <- names(env_layers())
-    }else{
+    reduce_colin <- names(env_layers())
+    updateSelectInput(inputId = "extract_variables", choices = names(env_layers()),
+                      selected = reduce_colin)
+
+    tryCatch({kept_colin_var$rc <- reduce_colin()[[3]]}, error = function(e)print(e))
+    if (!is.null(isolate(kept_colin_var$rc))) {
       reduce_colin <- names(env_layers())[!names(env_layers()) %in% reduce_colin()[[3]]]
+      updateSelectInput(inputId = "extract_variables", choices = names(env_layers()),
+                        selected = reduce_colin)
     }
 
-    updateSelectInput(inputId = "extract_variables", choices = names(env_layers()),
-                      selected = reduce_colin)
+
   })
   ## process extraction
   extracted_df <- eventReactive(input$extract_data, {
     req(valided_dp())
-    sdm_extract(
+    nm_extract(
       data = valided_dp(),
-      x = ifelse(input$partition_type %in% c("part_sband", "part_sblock", "part_senv"), "x", input$long_var),
-      y = ifelse(input$partition_type %in% c("part_sband", "part_sblock", "part_senv"), "y", input$lat_var),
+      longitude = ifelse(input$partition_type %in% c("part_sband", "part_sblock", "part_senv"), "x", input$long_var),
+      latitude = ifelse(input$partition_type %in% c("part_sband", "part_sblock", "part_senv"), "y", input$lat_var),
       env_layer = env_layers(),
+      current_crs = input$occ_current_crs,
       variables = input$extract_variables,
       filter_na = TRUE
     )
diff --git a/inst/nimo/ui.R b/inst/nimo/ui.R
index d785b1d..8a9036d 100644
--- a/inst/nimo/ui.R
+++ b/inst/nimo/ui.R
@@ -268,7 +268,10 @@ nimo_body <- shinydashboard::dashboardBody(
                                                   selected = "Presence"),
                                       shinyFilesButton("choose_data_file", "Load data",
                                                        "Select Species occurrence data", multiple = FALSE,
-                                                       icon = icon("file-upload"))
+                                                       icon = icon("file-upload")),
+                                      selectInput("occ_current_crs", label = "Coordinates CRS",
+                                                  choices = c("Geographic" = "geographic", "Projected" = "projected"), multiple = FALSE,
+                                                  selected = c("Geographic" = "geographic"))
                                       ),
                      hr(),
                      conditionalPanel("output.geo_distribution != null",
diff --git a/tests/testthat/test-nm_extract.R b/tests/testthat/test-nm_extract.R
new file mode 100644
index 0000000..c822630
--- /dev/null
+++ b/tests/testthat/test-nm_extract.R
@@ -0,0 +1,22 @@
+test_that("test for nm_extract", {
+  require(terra)
+
+  wbv_path <- paste0(system.file("extdata", package = "nimo"), "/WbV_subset_occ_KNP.csv")
+  wbv_df <- read.csv(wbv_path)
+
+  env_layers_path <- paste0(system.file("extdata", package = "nimo"), "/env_layers")
+  env_layers <- terra::rast(list.files(env_layers_path, full.names = TRUE))
+
+  extracted <- nimo::nm_extract(data = wbv_df,
+                         longitude = "decimalLongitude",
+                         latitude = "decimalLatitude",
+                         env_layer = env_layers)
+
+  expect_true(all(class(extracted) %in% c("tbl_df", "tbl", "data.frame")))
+  expect_error(nm_extract(data = wbv_df,
+                          longitude = "decimalLongitude",
+                          latitude = "decimalLatitude",
+                          env_layer = env_layers,
+                          current_crs = "geo"))
+  rm(extracted, wbv_path, wbv_df, env_layers_path, env_layers)
+})