diff --git a/DESCRIPTION b/DESCRIPTION
index 99fdd50..0f0b2a4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -32,7 +32,7 @@ Imports:
     mapdata,
     TMB,
     MatrixModels,
-    rgdal,
+    sf,
     ThorsonUtilities,
     abind,
     TMBhelper
@@ -47,6 +47,6 @@ Additional_repositories: https://www.math.ntnu.no/inla/R/stable
 License: GPL-2
 LazyData: yes
 BuildVignettes: yes
-RoxygenNote: 5.0.1
+RoxygenNote: 6.0.1
 URL: http://github.com/nwfsc-assess/geostatistical_delta-GLMM
 BugReports: http://github.com/nwfsc-assess/geostatistical_delta-GLMM/issues
diff --git a/R/Convert_LL_to_EastNorth_Fn.R b/R/Convert_LL_to_EastNorth_Fn.R
index 029224e..e682fd0 100644
--- a/R/Convert_LL_to_EastNorth_Fn.R
+++ b/R/Convert_LL_to_EastNorth_Fn.R
@@ -1,4 +1,3 @@
-
 #' Convert from Lat-Long to Eastings-Northings using WGS
 #'
 #' \code{Convert_LL_to_EastNorth_Fn} converts from Latitude-Longitude to World Geodetic System Eastings-Northings for a given location
@@ -6,7 +5,8 @@
 #' @param Lat vector of latitudes
 #' @param Lon vector of longitudes
 #' @param crs EPSG reference for coordinate reference system (CRS) defining Eastings-Northings after transformation
-
+#' @author James Thorson, Sophie Mormede
+#' 
 #' @return A data frame with the following columns
 #' \describe{
 #'   \item{E_km}{The eastings for each value of Lon (in kilometers)}
@@ -17,22 +17,25 @@
 Convert_LL_to_EastNorth_Fn <-
 function( Lon, Lat, crs=NA ){
   # SEE:  https://github.com/nwfsc-assess/geostatistical_delta-GLMM/issues/25#issuecomment-345825230
-
+  # SEE: https://github.com/nwfsc-assess/geostatistical_delta-GLMM/issues/38
+  
   # Attach package
-  require(rgdal)
-  on.exit( detach("package:rgdal") )
-
+  require(sf)
+  on.exit( detach("package:sf") )
+  
   # Transform
   dstart<-data.frame(lon=Lon, lat=Lat) # that's the object
-  coordinates(dstart) <- c("lon", "lat")
-  proj4string(dstart) <- CRS("+init=epsg:4326") # that's the lat long projection
-  CRS.new <- CRS(crs) # that's the eastings and northings projection
-  dstart.t <- spTransform(dstart, CRS.new) # here's where you transform
-
+  dstart<-as.matrix(dstart)
+  
+  dstart <- st_multipoint(dstart)
+  dstart <- st_sfc(dstart, crs = 4326)
+  dstart.t<- st_transform(dstart, crs)
+  dstart.t<-st_coordinates(dstart.t)
+  
   # Clean up
-  dstart.t = cbind( "E_km"=dstart.t@coords[,"lon"]/1000, "N_km"=dstart.t@coords[,'lat']/1000 )
+  dstart.t = cbind( "E_km"=dstart.t[,1]/1000, "N_km"=dstart.t[,2]/1000 )
   attr(dstart.t,"zone") = crs
-
+  
   # Return results
   return( dstart.t )
 }