diff --git a/DESCRIPTION b/DESCRIPTION
index d258984..1e91597 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: geoGraph
 Type: Package
 Title: Walking through the geographic space using graphs
-Version: 1.1.1.9001
+Version: 1.1.1.9002
 Authors@R: c(
     person("Thibaut", "Jombart", role = "aut"),
     person("Andrea", "Manica", email = "am315@cam.ac.uk", role = c("aut", "cre")))
@@ -23,14 +23,17 @@ Depends:
   methods,
   graph
 Imports:
-  fields,
-  RBGL,
-  sp
+    fields,
+    RBGL,
+    rnaturalearth,
+    rnaturalearthdata,
+    sp,
+    sf,
+    magrittr
 Suggests: 
   testthat, 
   knitr, 
-  rmarkdown,
-  sf
+  rmarkdown
 RoxygenNote: 7.2.3
 Collate: 
     'classes.R'
@@ -55,5 +58,6 @@ Collate:
     'rebuild.R'
     'setCosts.R'
     'setDistCosts.R'
+    'utils-pipe.R'
     'zoom.R'
     'zzz.R'
diff --git a/NAMESPACE b/NAMESPACE
index cdba679..ca52fd4 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(plot,gPath)
+export("%>%")
 export()
 export(.gData.valid)
 export(.gGraph.valid)
@@ -77,3 +78,4 @@ import(sp)
 importFrom(graphics,identify)
 importFrom(graphics,locator)
 importFrom(graphics,segments)
+importFrom(magrittr,"%>%")
diff --git a/R/classes.R b/R/classes.R
index 260b040..8e2ab5c 100644
--- a/R/classes.R
+++ b/R/classes.R
@@ -281,15 +281,43 @@ setMethod("initialize", "gGraph", function(.Object, ...) {
 
   ## handle @coords ##
   if (!is.null(input$coords)) {
-    if (is.list(input$coords) && length(input$coords) == 2) {
+
+
+    if (is.list(input$coords)) {
       input$coords <- as.data.frame(input$coords)
     }
 
     if (is.data.frame(input$coords)) {
       input$coords <- as.matrix(input$coords)
     }
-
-    if (nrow(input$coords) > 0 && !is.numeric(input$coords)) stop("Argument coords has to be numeric.")
+    
+    if (ncol(input$coords)!=2){
+      stop("Argument coords must include only two columns (longitude and latitude).")
+    }
+    
+    if (nrow(input$coords) > 0 && !is.numeric(input$coords)) {
+      stop("Argument coords has to be numeric.")
+      }
+    
+    ## NAs in coords
+    if (any(is.na(input$coords))) {
+      stop("Argument coords includes NAs")
+    }
+    
+    ## Convert all column names to lower case
+    colnames(input$coords) <- tolower(colnames(input$coords))
+    ## Create list of lon/lat column heading names
+    lonlist <- list("lon", "long", "longitude", "x")
+    latlist <- list("lat", "latitude", "y")
+    ## Test if the column order is inverted
+    if (is.element(colnames(input$coords)[1], latlist) & 
+        is.element(colnames(input$coords)[2], lonlist)) {
+      input$coords[, c(1, 2)] <- input$coords[, c(2, 1)]
+    } else if  (!(is.element(colnames(input$coords)[1], lonlist) & 
+                  is.element(colnames(input$coords)[2], latlist))){
+      message("The coordinate column names are not part of the standardised list;\n",
+      "we will use the order they were given in, make sure it corresponds to x and y!")
+    } # if neither of the if catches it, then the names are part of the lists and in the correct order
 
     ## names of the matrix
     colnames(input$coords) <- c("lon", "lat")
diff --git a/R/extractFromLayer.R b/R/extractFromLayer.R
index 7b4cb12..5ce28ae 100644
--- a/R/extractFromLayer.R
+++ b/R/extractFromLayer.R
@@ -49,15 +49,9 @@
 #'
 #' plot(worldgraph.10k, reset = TRUE)
 #'
-#'
-#' ## see what info is available
-#' names(worldshape@data)
-#' unique(worldshape@data$CONTINENT)
-#'
-#'
 #' ## retrieve continent info for all nodes
 #' ## (might take a few seconds)
-#' x <- extractFromLayer(worldgraph.10k, layer = worldshape, attr = "CONTINENT")
+#' x <- extractFromLayer(worldgraph.10k, layer = "world", attr = "CONTINENT")
 #' x
 #' table(getNodesAttr(x, attr.name = "CONTINENT"))
 #'
@@ -90,67 +84,61 @@ setGeneric("extractFromLayer", function(x, ...) {
 #' @rdname extractFromLayer
 #' @export
 setMethod("extractFromLayer", "matrix", function(x, layer = "world", attr = "all", ...) {
+  
+  
   ## This functions automatically assigns to land all points overlapping the country polygons
   # if(!require(maptools)) stop("maptools package is required.")
 
   ## Load default shapefile ##
   if (is.character(layer) && layer[1] == "world") {
-    layer <- worldshape
+    # use rnaturalearth instead of the inbuilt dataset
+    # layer <- rnaturalearth::ne_countries(scale="medium", returnclass = "sf")
+    # sf::sf_use_s2(FALSE)
+    layer <- sf::st_read(system.file("files/shapefiles/world-countries.shp", package = "geoGraph"))
   }
 
   ## TODO if the layer is null, we should throw an error!!!
   if (!is.null(layer)) {
-    if (!inherits(layer, "SpatialPolygonsDataFrame")) {
-      stop("Layer must be a SpatialPolygonsDataFrame object \n(see st_read and as_Spatial in sf to import such data from a GIS shapefile).")
+    if (!inherits(layer, "sf")) {
+      if (inherits(layer, "SpatialPolygonsDataFrame")){
+        layer <- sf::st_as_sf(layer)
+      } else {
+        stop("Layer must be a sf object \n(see st_read in sf to import such data from a GIS shapefile).")
+      }
     }
   }
 
+
   ## search attr in data ##
   if (attr[1] == "all") {
-    selAttr <- 1:ncol(layer@data)
+    selAttr <- 1:ncol(layer)
   } else {
-    selAttr <- match(attr, colnames(layer@data)) # selected attributes
+    selAttr <- match(attr, colnames(layer)) # selected attributes
     if (any(is.na(selAttr))) { # attribute not found in layer@data
       cat("\nSome requested attribute (attr) not found in the layer.\n")
       cat("\nAvailable data are:\n")
-      print(utils::head(layer@data))
+      print(utils::head(layer))
       return(NULL) # return NULL if attr not found, not generate an error
     }
   }
 
-  ## variables and initialization ##
-  long <- unlist(x[, 1]) # unlist needed when nrow==1
-  lat <- unlist(x[, 2])
-  n.poly.list <- length(layer@polygons) # number of lists of Polygons obj.
-  res <- NULL
-  dat <- layer@data
-  layerId <- rep(NA, length(long)) # stores the id of matching polygon for each location
-
 
-  ## main computations ##
+  # create an sf point object from the coordinates
+  locations_st <- x %>% as.data.frame %>% 
+    sf::st_as_sf(coords=c(1,2)) %>%
+    sf::st_set_crs(sf::st_crs(layer))
+  # now find points in polygons
+  points_within <- sf::st_intersects(layer, locations_st)
+  points_within <- data.frame(x = unlist(points_within), 
+                              polygon = rep(seq_along(lengths(points_within)), lengths(points_within)))
+  points_assignment <- data.frame(x=seq(1, nrow(x)), polygon = NA)
+  # add missing points for which we have no information
+  points_assignment[points_within$x,"polygon"]<-points_within$polygon
 
-  ## browsing elements of @polygons
-  ## each is a list with a @Polygons slot
-  for (i in 1:n.poly.list) {
-    this.poly.list <- layer@polygons[[i]]
-    n.polys <- length(this.poly.list@Polygons)
-    points.in.this.poly <- rep(0, length(long))
+  dat <- layer %>% sf::st_drop_geometry()
+  # @TOFIX the line below will fail if layerId is all NAs (i.e. no points were assigned to a polygon)
+  res <- dat[points_assignment$polygon, selAttr, drop = FALSE]
 
-    ## browsing elements of @Polygons
-    for (j in 1:n.polys) { ##
-      this.poly <- this.poly.list@Polygons[[j]]
-      points.in.this.poly <- points.in.this.poly +
-        sp::point.in.polygon(long, lat, this.poly@coords[, 1], this.poly@coords[, 2])
-
-      points.in.this.poly <- as.logical(points.in.this.poly)
-
-      if (any(points.in.this.poly)) {
-        layerId[points.in.this.poly] <- this.poly.list@ID
-      }
-    } # end for j
-  } # end for i
-
-  res <- dat[layerId, selAttr, drop = FALSE]
   row.names(res) <- rownames(x)
 
   return(res)
@@ -237,4 +225,4 @@ setMethod("extractFromLayer", "gData", function(x, layer = "world", attr = "all"
   }
 
   return(x)
-}) # end findLand
+})
diff --git a/R/findLand.R b/R/findLand.R
index 6de6300..bbbcca8 100644
--- a/R/findLand.R
+++ b/R/findLand.R
@@ -91,6 +91,10 @@ setMethod("findLand", "matrix", function(x, shape = "world", ...) {
     }
   }
 
+  if (any(is.na(x))) {
+    stop("Matrix contains NA values.")
+  }
+  
   long <- x[, 1]
   lat <- x[, 2]
   n.country <- length(shape@polygons)
diff --git a/R/utils-pipe.R b/R/utils-pipe.R
new file mode 100644
index 0000000..fd0b1d1
--- /dev/null
+++ b/R/utils-pipe.R
@@ -0,0 +1,14 @@
+#' Pipe operator
+#'
+#' See \code{magrittr::\link[magrittr:pipe]{\%>\%}} for details.
+#'
+#' @name %>%
+#' @rdname pipe
+#' @keywords internal
+#' @export
+#' @importFrom magrittr %>%
+#' @usage lhs \%>\% rhs
+#' @param lhs A value or the magrittr placeholder.
+#' @param rhs A function call using the magrittr semantics.
+#' @return The result of calling `rhs(lhs)`.
+NULL
diff --git a/man/extractFromLayer.Rd b/man/extractFromLayer.Rd
index 83b2f35..6947de3 100644
--- a/man/extractFromLayer.Rd
+++ b/man/extractFromLayer.Rd
@@ -71,15 +71,9 @@ input formats.
 
 plot(worldgraph.10k, reset = TRUE)
 
-
-## see what info is available
-names(worldshape@data)
-unique(worldshape@data$CONTINENT)
-
-
 ## retrieve continent info for all nodes
 ## (might take a few seconds)
-x <- extractFromLayer(worldgraph.10k, layer = worldshape, attr = "CONTINENT")
+x <- extractFromLayer(worldgraph.10k, layer = "world", attr = "CONTINENT")
 x
 table(getNodesAttr(x, attr.name = "CONTINENT"))
 
diff --git a/man/pipe.Rd b/man/pipe.Rd
new file mode 100644
index 0000000..a648c29
--- /dev/null
+++ b/man/pipe.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils-pipe.R
+\name{\%>\%}
+\alias{\%>\%}
+\title{Pipe operator}
+\usage{
+lhs \%>\% rhs
+}
+\arguments{
+\item{lhs}{A value or the magrittr placeholder.}
+
+\item{rhs}{A function call using the magrittr semantics.}
+}
+\value{
+The result of calling \code{rhs(lhs)}.
+}
+\description{
+See \code{magrittr::\link[magrittr:pipe]{\%>\%}} for details.
+}
+\keyword{internal}
diff --git a/tests/testthat/test_classes.R b/tests/testthat/test_classes.R
index e019f46..986a855 100644
--- a/tests/testthat/test_classes.R
+++ b/tests/testthat/test_classes.R
@@ -1,10 +1,70 @@
-library("geoGraph")
-context("Test classes")
 
 test_that("Empty constructors work", {
-  skip_on_cran()
   x <- new("gGraph")
+  expect_true(inherits(x,"gGraph"))
   y <- new("gData")
-  expect_is(x, "gGraph")
-  expect_is(y, "gData")
+  expect_true(inherits(y,"gData"))
 })
+
+test_that("Contructors fails with invalid coordinates",{
+  NACoords <- data.frame(long = c(-24, NA), lat =  c(31,55))
+  #Create gGraph with NA's
+  expect_error(new("gGraph", coords = NACoords),
+               "Argument coords includes NAs")
+
+})
+
+test_that("Contructors fails with invalid matrix dimensions",{
+  extra_coords <- data.frame(long = c(-24, 37), lat =  c(31,55), x = c(31,31))
+  #Create gGraph with three columns
+  expect_error(new("gGraph", coords = extra_coords),
+               "Argument coords must include")
+
+})
+
+test_that("Contructors fails with invalid non-numeric matrix",{
+  non_num_coords <- data.frame(long = c("lon1", 37), lat =  c(31,55))
+  #Create gGraph with non numeric elements
+  expect_error(new("gGraph", coords = non_num_coords),
+               "Argument coords has to be numeric")
+  
+})
+
+
+
+test_that("Constructor accounts for different column names" , {
+  columnheading_names <-
+    data.frame(longitude =  c(31, 55), Latitude = c(-24, 37))
+  columnheading_names <- new("gGraph", coords = columnheading_names)
+  new_columnheading_names <-
+    data.frame(lon = c(31, 55), lat =  c(-24, 37))
+  new_columnheading_names <-
+    new("gGraph", coords = new_columnheading_names)
+  expect_identical(columnheading_names,
+                   new_columnheading_names)
+  
+})
+
+test_that("Constructor reverses coord column order" , {
+  column_heading <- data.frame(lat = c(-24, 37), lon =  c(31,55)) 
+  #Create Ggraph with lat/lon headings
+  correct_heading <- new("gGraph", coords = column_heading)
+  column_heading <- data.frame(lon =  c(31,55), lat = c(-24, 37))
+  #Create Ggraph with lon/lat headings
+  swapped_heading <- new("gGraph", coords = column_heading)
+  expect_identical(correct_heading, swapped_heading)
+})
+
+test_that("we give message when columns are not recognised",{
+  column_heading <- data.frame(lon=  c(31,55), lat = c(-24, 37))
+  #Create Ggraph with lat/lon headings
+  correct_heading <- new("gGraph", coords = column_heading)
+  column_heading <- data.frame(blah=  c(31,55), lat = c(-24, 37))  
+  #Create Ggraph with lon/lat headings
+  expect_message(unrecognised_heading <- new("gGraph", coords = column_heading),
+                 "The coordinate column names are not part of the standardised list")
+  expect_identical(correct_heading, unrecognised_heading)
+}
+          
+          
+          )
diff --git a/tests/testthat/test_extractFromLayer.R b/tests/testthat/test_extractFromLayer.R
new file mode 100644
index 0000000..b0cad72
--- /dev/null
+++ b/tests/testthat/test_extractFromLayer.R
@@ -0,0 +1,8 @@
+test_that("extractFromLayer assigns points correctly",
+          {
+            # create a matrix of locations including two continents and the sea
+            myCoords <- data.frame(long = c(-24, 71.5, -46.5), lat =  c(31, 30,-23.5))
+            # assign to continents
+            continents <- extractFromLayer(myCoords, layer = worldshape, attr = "CONTINENT")
+            expect_identical(as.character(continents$CONTINENT), c(NA,"Asia","South America"))
+          })
\ No newline at end of file
diff --git a/tests/testthat/test_findLand.R b/tests/testthat/test_findLand.R
new file mode 100644
index 0000000..17af0be
--- /dev/null
+++ b/tests/testthat/test_findLand.R
@@ -0,0 +1,48 @@
+test_that("find land correctly",{
+  # create a gGraph with one sea node and one land node
+  myCoords <- data.frame(lon = c(-24, 37), lat =  c(31,55))
+  obj <- new("gGraph", coords = myCoords)
+  obj <- findLand(obj) 
+  # check that we classified the nodes correctly
+  expect_true(all(obj@nodes.attr$habitat==c("sea","land")))
+  
+  # error if we pass an incorrect class
+  expect_error(findLand("blah"),
+               "unable to find an inherited method")
+  
+  # tests error for NAs as coordinates (this only applies to matrices and data.frames)
+  # gGraph objects should fail when creted with NA coordinates
+  NACoords <- data.frame(long = c(-24, NA), lat =  c(31,55))
+  expect_error(findLand(NACoords),
+               "Matrix contains NA")
+  
+ 
+
+})
+
+
+#Test whether findLand generates an error with invalid matrix object
+test_that("co-ordinate format",{
+  #Create co-ordinates matrix
+  obj <- matrix(c(-24, 31, 37, 55), nrow=2,ncol=2, byrow=TRUE)
+  obj <- findLand(obj) 
+  
+  #check for factor output
+  expect_is(obj,"factor")
+  #check for correct output
+  expect_equal(obj,factor(c('sea','land')))
+  
+  #Create co-ordinates matrix with NA
+  NA_matrix <- matrix(c(-24, NA, 37, 55), nrow=2,ncol=2, byrow=TRUE)
+  #NA entries are recognised and produce error in plot 
+  expect_true(is.na(NA_matrix[1,2]))
+  expect_error(plotEdges(NA_matrix))
+  
+  #NA produces error in findLand
+  expect_error(findLand(NA_matrix))
+  
+  
+})
+
+
+