diff --git a/R/cc_outl.R b/R/cc_outl.R
index 93b57fe..9f15756 100644
--- a/R/cc_outl.R
+++ b/R/cc_outl.R
@@ -199,6 +199,13 @@ cc_outl <- function(x,
                            area = terra::expanse(ref))
         area <- area[!is.na(area$area), ]
         area <- area[!is.na(area$country), ]
+
+        ## naturalearth uses -99 for missing values?
+        area <- area[area$country != "-99", ]
+
+        ## CN-TW (Taiwan) doesn't map to a country in the GBIF database:
+        ## naturalearth dataset uses CN-TW for Taiwan, GBIF uses TW.
+        area[area$country == "CN-TW", "country"] <- "TW"
         
         # get number of records in GBIF per country as proxy for sampling intensity
         nrec <- vapply(area$country, 
diff --git a/R/cc_sea.R b/R/cc_sea.R
index 6a3dbc2..d52cc2c 100644
--- a/R/cc_sea.R
+++ b/R/cc_sea.R
@@ -137,7 +137,9 @@ cc_sea <- function(x,
     ## have to make sure they are the same and in this case they already are
     ## -----
     ## point-in-polygon test
-    ext_dat <- terra::extract(ref, pts)
+    ##ext_dat <- terra::extract(ref, pts)
+    ext_dat <- terra::relate(pts, ref, "intersects", pairs = TRUE,
+                             na.rm = FALSE)
     out <- !is.na(ext_dat[!duplicated(ext_dat[, 1]), 2])
     out <- data.frame(terra::geom(pts)[, c("x", "y"), drop = FALSE], out)
     colnames(out)[1:2] <- c(lon, lat)
@@ -154,7 +156,9 @@ cc_sea <- function(x,
     
     # select relevant columns
     #point in polygon test
-    ext_dat <- terra::extract(ref, pts)
+    ext_dat <- terra::relate(pts, ref, "intersects", pairs = TRUE,
+                             na.rm = FALSE) 
+    ##ext_dat <- terra::extract(ref, pts)
     out <- !is.na(ext_dat[!duplicated(ext_dat[, 1]), 2])
   }