diff --git a/DESCRIPTION b/DESCRIPTION
index e11ae10f..537ce3a7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: spdep
 Version: 1.3-5
-Date: 2024-06-02
+Date: 2024-06-10
 Title: Spatial Dependence: Weighting Schemes, Statistics
 Encoding: UTF-8
 Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"), 
@@ -38,7 +38,7 @@ Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"),
 	person("René", "Westerholt", role="ctb"),
 	person("Levi", "Wolf", role = "ctb"),
 	person("Danlin", "Yu", role = "ctb"))
-Depends: R (>= 3.3.0), methods, spData (>= 0.2.6.0), sf 
+Depends: R (>= 3.3.0), methods, spData (>= 2.3.1), sf 
 Imports: stats, deldir, boot (>= 1.3-1), graphics, utils, grDevices, units, s2, e1071, sp (>= 1.0)
 Suggests: spatialreg (>= 1.2-1), Matrix, parallel, dbscan, RColorBrewer, lattice, xtable, foreign, igraph, RSpectra, knitr, classInt, tmap, spam, ggplot2, rmarkdown, tinytest, rgeoda
 URL: https://github.com/r-spatial/spdep/, https://r-spatial.github.io/spdep/
diff --git a/NEWS.md b/NEWS.md
index b079dc28..8521f627 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,6 @@
-# Version 1.3-4 (development)
+# Version 1.3-5 (development)
+
+* #157 migrate ESRI Shapefile to GPKG files; convert bhicv.shp to GPKG 
 
 * #155 Throw error if `hotspot` despatched on object without a `"quadr"` attribute
 
diff --git a/inst/etc/shapes/bhicv.dbf b/inst/etc/shapes/bhicv.dbf
deleted file mode 100644
index a1de0076..00000000
Binary files a/inst/etc/shapes/bhicv.dbf and /dev/null differ
diff --git a/inst/etc/shapes/bhicv.gpkg b/inst/etc/shapes/bhicv.gpkg
new file mode 100644
index 00000000..429d0685
Binary files /dev/null and b/inst/etc/shapes/bhicv.gpkg differ
diff --git a/inst/etc/shapes/bhicv.shp b/inst/etc/shapes/bhicv.shp
deleted file mode 100644
index 927d1882..00000000
Binary files a/inst/etc/shapes/bhicv.shp and /dev/null differ
diff --git a/inst/etc/shapes/bhicv.shx b/inst/etc/shapes/bhicv.shx
deleted file mode 100644
index ffc19647..00000000
Binary files a/inst/etc/shapes/bhicv.shx and /dev/null differ
diff --git a/inst/etc/shapes/columbus.dbf b/inst/etc/shapes/columbus.dbf
deleted file mode 100644
index f850bda2..00000000
Binary files a/inst/etc/shapes/columbus.dbf and /dev/null differ
diff --git a/inst/etc/shapes/columbus.shp b/inst/etc/shapes/columbus.shp
deleted file mode 100644
index d6dde43c..00000000
Binary files a/inst/etc/shapes/columbus.shp and /dev/null differ
diff --git a/inst/etc/shapes/columbus.shx b/inst/etc/shapes/columbus.shx
deleted file mode 100644
index 7086b6f2..00000000
Binary files a/inst/etc/shapes/columbus.shx and /dev/null differ
diff --git a/inst/etc/shapes/eire.dbf b/inst/etc/shapes/eire.dbf
deleted file mode 100644
index 14f2d8e1..00000000
Binary files a/inst/etc/shapes/eire.dbf and /dev/null differ
diff --git a/inst/etc/shapes/eire.shp b/inst/etc/shapes/eire.shp
deleted file mode 100644
index 612397bd..00000000
Binary files a/inst/etc/shapes/eire.shp and /dev/null differ
diff --git a/inst/etc/shapes/eire.shx b/inst/etc/shapes/eire.shx
deleted file mode 100644
index 580107d7..00000000
Binary files a/inst/etc/shapes/eire.shx and /dev/null differ
diff --git a/man/EBImoran.mc.Rd b/man/EBImoran.mc.Rd
index 88fee755..b6699526 100644
--- a/man/EBImoran.mc.Rd
+++ b/man/EBImoran.mc.Rd
@@ -60,7 +60,7 @@ following components:
 \seealso{\code{\link{moran}}, \code{\link{moran.mc}}, \code{\link{EBest}}}
 
 \examples{
-nc.sids <- st_read(system.file("shapes/sids.shp", package="spData")[1], quiet=TRUE)
+nc.sids <- st_read(system.file("shapes/sids.gpkg", package="spData")[1], quiet=TRUE)
 rn <- as.character(nc.sids$FIPS)
 ncCC89_nb <- read.gal(system.file("weights/ncCC89.gal", package="spData")[1],
  region.id=rn)
diff --git a/man/EBest.Rd b/man/EBest.Rd
index 98a86175..12c222df 100644
--- a/man/EBest.Rd
+++ b/man/EBest.Rd
@@ -32,7 +32,7 @@ and a \code{parameters} attribute list with components:
 \seealso{\code{\link{EBlocal}}, \code{\link{probmap}}, \code{\link{EBImoran.mc}}}
 
 \examples{
-auckland <- st_read(system.file("shapes/auckland.shp", package="spData")[1], quiet=TRUE)
+auckland <- st_read(system.file("shapes/auckland.gpkg", package="spData")[1], quiet=TRUE)
 res <- EBest(auckland$M77_85, 9*auckland$Und5_81)
 attr(res, "parameters")
 auckland$estmm000 <- res$estmm*1000
diff --git a/man/EBlocal.Rd b/man/EBlocal.Rd
index f08d92e4..b953add6 100644
--- a/man/EBlocal.Rd
+++ b/man/EBlocal.Rd
@@ -37,7 +37,7 @@ and a \code{parameters} attribute list with components (if both are zero, the es
 \seealso{\code{\link{EBest}}, \code{\link{probmap}}}
 
 \examples{
-auckland <- st_read(system.file("shapes/auckland.shp", package="spData")[1], quiet=TRUE)
+auckland <- st_read(system.file("shapes/auckland.gpkg", package="spData")[1], quiet=TRUE)
 auckland.nb <- poly2nb(auckland)
 res <- EBlocal(auckland$M77_85,  9*auckland$Und5_81, auckland.nb)
 auckland$est000 <- res$est*1000
diff --git a/man/SD.RStests.Rd b/man/SD.RStests.Rd
index f5260b17..e6aa317f 100644
--- a/man/SD.RStests.Rd
+++ b/man/SD.RStests.Rd
@@ -42,7 +42,7 @@ The results in the example below agree with those in Table 3, p. 22 in Koley and
 \code{\link{lm}}, \code{\link{lm.RStests}}
 }
 \examples{
-columbus <- sf::st_read(system.file("shapes/columbus.shp", package="spData")[1])
+columbus <- sf::st_read(system.file("shapes/columbus.gpkg", package="spData")[1])
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 col.listw <- nb2listw(col.gal.nb, style="W")
 lm_obj <- lm(CRIME ~ INC + HOVAL, data=columbus)
diff --git a/man/autocov_dist.Rd b/man/autocov_dist.Rd
index 49346e1c..9017a8d1 100644
--- a/man/autocov_dist.Rd
+++ b/man/autocov_dist.Rd
@@ -35,7 +35,7 @@ for details.}
 
 \seealso{\code{\link{nb2listw}}}
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 #xy <- cbind(columbus$X, columbus$Y)
 xy <- st_coordinates(st_centroid(st_geometry(columbus),
  of_largest_polygon=TRUE))
diff --git a/man/bhicv.Rd b/man/bhicv.Rd
index c6117fff..fd9b3a35 100644
--- a/man/bhicv.Rd
+++ b/man/bhicv.Rd
@@ -24,7 +24,7 @@
 %%\format{}
 %%\details{}
 \examples{
-bh <- st_read(system.file("etc/shapes/bhicv.shp",
+bh <- st_read(system.file("etc/shapes/bhicv.gpkg",
       package="spdep")[1])
 }
 \keyword{data}% at least one, from doc/KEYWORDS
diff --git a/man/choynowski.Rd b/man/choynowski.Rd
index 7fb00bab..fe063a4f 100644
--- a/man/choynowski.Rd
+++ b/man/choynowski.Rd
@@ -31,7 +31,7 @@ Interactive Spatial Data Analysis, Harlow: Longman, pp. 300--303.}
 
 \seealso{\code{\link{probmap}}}
 \examples{
-auckland <- st_read(system.file("shapes/auckland.shp", package="spData")[1], quiet=TRUE)
+auckland <- st_read(system.file("shapes/auckland.gpkg", package="spData")[1], quiet=TRUE)
 auckland.nb <- poly2nb(auckland)
 res <- choynowski(auckland$M77_85, 9*auckland$Und5_81)
 resl <- choynowski(auckland$M77_85, 9*auckland$Und5_81, legacy=TRUE)
diff --git a/man/columbus.Rd b/man/columbus.Rd
index 3e43c3ad..69d65869 100644
--- a/man/columbus.Rd
+++ b/man/columbus.Rd
@@ -13,7 +13,7 @@
 The data set is now part of the spData package
 }
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 }
 
diff --git a/man/compon.Rd b/man/compon.Rd
index 9e505e00..fe6c62d5 100644
--- a/man/compon.Rd
+++ b/man/compon.Rd
@@ -22,7 +22,7 @@ A list of:
 \author{Nicholas Lewin-Koh \email{nikko@hailmail.net}}
 \seealso{\code{\link{plot.nb}}}
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 coords <- st_coordinates(st_centroid(st_geometry(columbus)))
 plot(col.gal.nb, coords, col="grey")
diff --git a/man/diffnb.Rd b/man/diffnb.Rd
index 3fbe8ba2..5d27e418 100644
--- a/man/diffnb.Rd
+++ b/man/diffnb.Rd
@@ -21,7 +21,7 @@ A neighbours list with class \code{nb}
 \author{Roger Bivand \email{Roger.Bivand@nhh.no}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 coords <- st_centroid(st_geometry(columbus), of_largest_polygon=TRUE)
 rn <- row.names(columbus)
 knn1 <- knearneigh(coords, 1)
diff --git a/man/dnearneigh.Rd b/man/dnearneigh.Rd
index f6452e9f..59ab93a2 100644
--- a/man/dnearneigh.Rd
+++ b/man/dnearneigh.Rd
@@ -32,7 +32,7 @@ for neighbours satisfying the distance criteria. See \code{\link{card}} for deta
 \seealso{\code{\link{knearneigh}}, \code{\link{card}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 coords <- st_centroid(st_geometry(columbus), of_largest_polygon=TRUE)
 rn <- row.names(columbus)
 k1 <- knn2nb(knearneigh(coords))
diff --git a/man/edit.nb.Rd b/man/edit.nb.Rd
index ba7817df..ebd8dcd3 100644
--- a/man/edit.nb.Rd
+++ b/man/edit.nb.Rd
@@ -31,7 +31,7 @@ The function returns an object of class \code{nb} with the edited list of intege
 
 \examples{
 \dontrun{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 class(columbus)
 if (FALSE) nnb1 <- edit.nb(col.gal.nb, polys=as(columbus, "Spatial"))
 }
diff --git a/man/globalG.test.Rd b/man/globalG.test.Rd
index 67cc6f5d..55fb932e 100644
--- a/man/globalG.test.Rd
+++ b/man/globalG.test.Rd
@@ -42,7 +42,7 @@ p. 276; Bivand RS, Wong DWS 2018 Comparing implementations of global and local i
 \seealso{\code{\link{localG}}}
 
 \examples{
-nc.sids <- st_read(system.file("shapes/sids.shp", package="spData")[1], quiet=TRUE)
+nc.sids <- st_read(system.file("shapes/sids.gpkg", package="spData")[1], quiet=TRUE)
 sidsrate79 <- (1000*nc.sids$SID79)/nc.sids$BIR79
 dists <- c(10, 20, 30, 33, 40, 50, 60, 70, 80, 90, 100)
 ndists <- length(dists)
diff --git a/man/graphneigh.Rd b/man/graphneigh.Rd
index 499c1126..9f9bc35d 100644
--- a/man/graphneigh.Rd
+++ b/man/graphneigh.Rd
@@ -103,7 +103,7 @@ The helper functions return an \code{nb} object with a list of integer
 \code{\link{knn2nb}}, \code{\link{card}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 sf_obj <- st_centroid(st_geometry(columbus), of_largest_polygon)
 sp_obj <- as(sf_obj, "Spatial")
 coords <- st_coordinates(sf_obj)
diff --git a/man/include.self.Rd b/man/include.self.Rd
index b42c40f7..465c784f 100644
--- a/man/include.self.Rd
+++ b/man/include.self.Rd
@@ -21,7 +21,7 @@ The function returns an object of class \code{nb} with a list of integer vectors
 \seealso{\code{\link{summary.nb}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 coords <- st_coordinates(st_centroid(columbus))
 summary(col.gal.nb, coords)
diff --git a/man/joincount.multi.Rd b/man/joincount.multi.Rd
index a2674389..08f8b187 100644
--- a/man/joincount.multi.Rd
+++ b/man/joincount.multi.Rd
@@ -34,7 +34,7 @@ A matrix with class \code{jcmulti} with row and column names for observed and ex
 \seealso{\code{\link{joincount.test}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 HICRIME <- cut(columbus$CRIME, breaks=c(0,35,80), labels=c("low","high"))
 (nb <- poly2nb(columbus))
 lw <- nb2listw(nb, style="B")
diff --git a/man/joincount.test.Rd b/man/joincount.test.Rd
index 6c8f77be..77f013c8 100644
--- a/man/joincount.test.Rd
+++ b/man/joincount.test.Rd
@@ -62,7 +62,7 @@ print(is.symmetric.nb(COL.k4.nb))
 joincount.test(HICRIME, nb2listw(COL.k4.nb, style="B"))
 cat("Note non-symmetric weights matrix - use listw2U()\n")
 joincount.test(HICRIME, listw2U(nb2listw(COL.k4.nb, style="B")))
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 HICRIME <- cut(columbus$CRIME, breaks=c(0,35,80), labels=c("low","high"))
 (nb <- poly2nb(columbus))
 lw <- nb2listw(nb, style="B")
diff --git a/man/knearneigh.Rd b/man/knearneigh.Rd
index cc968716..4efd43d3 100644
--- a/man/knearneigh.Rd
+++ b/man/knearneigh.Rd
@@ -32,7 +32,7 @@ A list of class \code{knn}
 \code{\link{knn2nb}}, \code{\link[dbscan]{kNN}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 coords <- st_centroid(st_geometry(columbus), of_largest_polygon=TRUE)
 col.knn <- knearneigh(coords, k=4)
 plot(st_geometry(columbus), border="grey")
diff --git a/man/knn2nb.Rd b/man/knn2nb.Rd
index b4ec8226..b4d87967 100644
--- a/man/knn2nb.Rd
+++ b/man/knn2nb.Rd
@@ -24,7 +24,7 @@ The function returns an object of class \code{nb} with a list of integer vectors
 \seealso{\code{\link{knearneigh}}, \code{\link{card}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 coords <- st_coordinates(st_centroid(columbus))
 col.knn <- knearneigh(coords, k=4)
 plot(st_geometry(columbus), border="grey")
diff --git a/man/listw2sn.Rd b/man/listw2sn.Rd
index a237ca5c..74dcef21 100644
--- a/man/listw2sn.Rd
+++ b/man/listw2sn.Rd
@@ -34,7 +34,7 @@ sn2listw(sn, style = NULL, zero.policy = NULL, from_mat2listw=FALSE)
 \seealso{\code{\link{nb2listw}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 col.listw <- nb2listw(col.gal.nb)
 col.listw$neighbours[[1]]
diff --git a/man/lm.morantest.exact.Rd b/man/lm.morantest.exact.Rd
index 060cb2a1..d0c744fa 100644
--- a/man/lm.morantest.exact.Rd
+++ b/man/lm.morantest.exact.Rd
@@ -52,7 +52,7 @@ must be one of greater (default), less or two.sided.}
 
 \seealso{\code{\link{lm.morantest.sad}}}
 \examples{
-eire <- st_read(system.file("shapes/eire.shp", package="spData")[1])
+eire <- st_read(system.file("shapes/eire.gpkg", package="spData")[1])
 row.names(eire) <- as.character(eire$names)
 st_crs(eire) <- "+proj=utm +zone=30 +ellps=airy +units=km"
 eire.nb <- poly2nb(eire)
diff --git a/man/lm.morantest.sad.Rd b/man/lm.morantest.sad.Rd
index 642be0cd..5aad56a7 100644
--- a/man/lm.morantest.sad.Rd
+++ b/man/lm.morantest.sad.Rd
@@ -68,7 +68,7 @@ Geographical Analysis, 34, pp. 187--206; Bivand RS, Wong DWS 2018 Comparing impl
 \seealso{\code{\link{lm.morantest}}}
 
 \examples{
-eire <- st_read(system.file("shapes/eire.shp", package="spData")[1])
+eire <- st_read(system.file("shapes/eire.gpkg", package="spData")[1])
 row.names(eire) <- as.character(eire$names)
 st_crs(eire) <- "+proj=utm +zone=30 +ellps=airy +units=km"
 eire.nb <- poly2nb(eire)
diff --git a/man/localGS.Rd b/man/localGS.Rd
index 8252e639..087897eb 100644
--- a/man/localGS.Rd
+++ b/man/localGS.Rd
@@ -38,7 +38,7 @@ Sonea, A. and Westerholt, R. (2021): Geographic and temporal access to basic ban
   
 \examples{
 \donttest{
-boston.tr <- sf::st_read(system.file("shapes/boston_tracts.shp", package="spData")[1])
+boston.tr <- sf::st_read(system.file("shapes/boston_tracts.gpkg", package="spData")[1])
 boston.tr_utm <- st_transform(boston.tr, 32619) #26786
 
 boston_listw1 <- nb2listwdist(dnearneigh(st_centroid(boston.tr_utm), 1, 2000),
diff --git a/man/localmoran.exact.Rd b/man/localmoran.exact.Rd
index 56d5d222..ecc2fc6f 100644
--- a/man/localmoran.exact.Rd
+++ b/man/localmoran.exact.Rd
@@ -72,7 +72,7 @@ When the alternative is being tested, a list of left and right M products in att
 
 \seealso{\code{\link{lm.morantest.exact}}, \code{\link{localmoran.sad}}}
 \examples{
-eire <- st_read(system.file("shapes/eire.shp", package="spData")[1])
+eire <- st_read(system.file("shapes/eire.gpkg", package="spData")[1])
 row.names(eire) <- as.character(eire$names)
 st_crs(eire) <- "+proj=utm +zone=30 +ellps=airy +units=km"
 eire.nb <- poly2nb(eire)
diff --git a/man/localmoran.sad.Rd b/man/localmoran.sad.Rd
index e6a7f607..d74b486b 100644
--- a/man/localmoran.sad.Rd
+++ b/man/localmoran.sad.Rd
@@ -105,7 +105,7 @@ Geographical Analysis, 34, pp. 187--206.}
  \code{\link{lm.morantest.sad}}, \code{\link{errorsarlm}}}
 
 \examples{
-eire <- st_read(system.file("shapes/eire.shp", package="spData")[1])
+eire <- st_read(system.file("shapes/eire.gpkg", package="spData")[1])
 row.names(eire) <- as.character(eire$names)
 st_crs(eire) <- "+proj=utm +zone=30 +ellps=airy +units=km"
 eire.nb <- poly2nb(eire)
diff --git a/man/localmoran_bv.Rd b/man/localmoran_bv.Rd
index eeba789b..52adfb0c 100644
--- a/man/localmoran_bv.Rd
+++ b/man/localmoran_bv.Rd
@@ -37,7 +37,7 @@ I_i^B= cx_i\Sigma_j{w_{ij}y_j}
 }
 \examples{
 # load columbus data
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData"))
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData"))
 nb <- poly2nb(columbus)
 listw <- nb2listw(nb)
 set.seed(1)
diff --git a/man/mat2listw.Rd b/man/mat2listw.Rd
index 9aad3c58..b686d48f 100644
--- a/man/mat2listw.Rd
+++ b/man/mat2listw.Rd
@@ -33,7 +33,7 @@ A \code{listw} object with the following members:
  
 \seealso{\code{\link{nb2listw}}, \code{\link{nb2mat}}}
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col005 <- dnearneigh(st_coordinates(st_centroid(st_geometry(columbus),
  of_largest_polygon=TRUE)), 0, 0.5, as.character(columbus$NEIGNO))
 summary(col005)
diff --git a/man/moran.test.Rd b/man/moran.test.Rd
index 7b57a3c4..e4017f20 100644
--- a/man/moran.test.Rd
+++ b/man/moran.test.Rd
@@ -83,7 +83,7 @@ moran.test(crime, nb2listw(COL.nb, style="W"), zero.policy=TRUE,
 moran.test(crime, nb2listw(COL.nb, style="W"), zero.policy=TRUE,
  na.action=na.exclude)
 moran.test(crime, nb2listw(COL.nb, style="W"), na.action=na.pass)
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col_geoms <- st_geometry(columbus)
 col_geoms[1] <- st_buffer(col_geoms[1], dist=-0.05)
 st_geometry(columbus) <- col_geoms
diff --git a/man/mstree.Rd b/man/mstree.Rd
index 81e2cb7a..c1b57504 100644
--- a/man/mstree.Rd
+++ b/man/mstree.Rd
@@ -45,7 +45,7 @@ mstree(nbw, ini = NULL)
 %%\seealso{ ~~objects to See Also as \code{\link{help}}, ~~~ }
 \examples{
 ### loading data
-bh <- st_read(system.file("etc/shapes/bhicv.shp",
+bh <- st_read(system.file("etc/shapes/bhicv.gpkg",
       package="spdep")[1], quiet=TRUE)
 st_crs(bh) <- "OGC:CRS84"
 ### data padronized
diff --git a/man/nb2lines.Rd b/man/nb2lines.Rd
index 843c7050..4268cc21 100644
--- a/man/nb2lines.Rd
+++ b/man/nb2lines.Rd
@@ -40,7 +40,7 @@ df2sn(df, i="i", i_ID="i_ID", j="j", wt="wt")
 \seealso{\code{\link{sn2listw}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 res <- listw2lines(nb2listw(col.gal.nb), st_geometry(columbus))
 summary(res)
diff --git a/man/nb2listw.Rd b/man/nb2listw.Rd
index 197e3264..a0e4bd88 100644
--- a/man/nb2listw.Rd
+++ b/man/nb2listw.Rd
@@ -42,7 +42,7 @@ and attributes:
 \seealso{\code{\link{summary.nb}}, \code{\link{read.gal}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 coords <- st_coordinates(st_centroid(columbus))
 cards <- card(col.gal.nb)
diff --git a/man/nb2mat.Rd b/man/nb2mat.Rd
index ed449af3..7cd71ef3 100644
--- a/man/nb2mat.Rd
+++ b/man/nb2mat.Rd
@@ -31,7 +31,7 @@ An n by n matrix, where n=length(neighbours)
 \seealso{\code{\link{nb2listw}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col005 <- dnearneigh(st_coordinates(st_centroid(st_geometry(columbus),
  of_largest_polygon=TRUE)), 0, 0.5, as.character(columbus$NEIGNO))
 summary(col005)
diff --git a/man/nbdists.Rd b/man/nbdists.Rd
index a5ac3884..c93caa37 100644
--- a/man/nbdists.Rd
+++ b/man/nbdists.Rd
@@ -24,7 +24,7 @@ A list with class \code{nbdist}
 \seealso{\code{\link{summary.nb}}, \code{\link{nb2listw}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 coords <- st_coordinates(st_centroid(columbus))
 dlist <- nbdists(col.gal.nb, coords)
diff --git a/man/nblag.Rd b/man/nblag.Rd
index 3dd7641a..215f5c41 100644
--- a/man/nblag.Rd
+++ b/man/nblag.Rd
@@ -23,7 +23,7 @@ returns a list of lagged neighbours lists each with class \code{nb}
 \seealso{\code{\link{summary.nb}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 coords <- st_coordinates(st_centroid(columbus))
 summary(col.gal.nb, coords)
diff --git a/man/nboperations.Rd b/man/nboperations.Rd
index a8f9a39b..0998a679 100644
--- a/man/nboperations.Rd
+++ b/man/nboperations.Rd
@@ -39,7 +39,7 @@ complement.nb(nb.obj)
 \code{\link{setdiff.nb}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 coords <- st_coordinates(st_centroid(columbus))
 col.tri.nb <- tri2nb(coords)
diff --git a/man/plot.nb.Rd b/man/plot.nb.Rd
index 9a9168df..442dca16 100644
--- a/man/plot.nb.Rd
+++ b/man/plot.nb.Rd
@@ -30,7 +30,7 @@ neighbours lists}
 \seealso{\code{\link{summary.nb}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 plot(col.gal.nb, st_geometry(columbus))
 title(main="GAL order 1 links with first nearest neighbours in red", cex.main=0.6)
diff --git a/man/poly2nb.Rd b/man/poly2nb.Rd
index 7230b377..8e45adb9 100644
--- a/man/poly2nb.Rd
+++ b/man/poly2nb.Rd
@@ -32,7 +32,7 @@ From 1.1-9 with \pkg{sf} 1.0-0, \pkg{s2} is used in bounding box indexing intern
 \seealso{\code{\link{summary.nb}}, \code{\link{card}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 coords <- st_coordinates(st_centroid(st_geometry(columbus)))
 xx <- poly2nb(as(columbus, "Spatial"))
@@ -62,7 +62,7 @@ fg[maxconts] <- "red"
 fg[xx[[maxconts]]] <- "green"
 plot(st_geometry(columbus), col=fg)
 title(main="Region with largest number of contiguities", cex.main=0.6)
-nc.sids <- st_read(system.file("shapes/sids.shp", package="spData")[1], quiet=TRUE)
+nc.sids <- st_read(system.file("shapes/sids.gpkg", package="spData")[1], quiet=TRUE)
 system.time(xxnb <- poly2nb(nc.sids))
 system.time(xxnb <- poly2nb(as(nc.sids, "Spatial")))
 plot(st_geometry(nc.sids))
diff --git a/man/probmap.Rd b/man/probmap.Rd
index 0afb68c1..c7ade839 100644
--- a/man/probmap.Rd
+++ b/man/probmap.Rd
@@ -29,7 +29,7 @@ probmap(n, x, row.names=NULL, alternative="less")
 \seealso{\code{\link{EBest}}, \code{\link{EBlocal}}, \code{\link{ppois}}}
 
 \examples{
-auckland <- st_read(system.file("shapes/auckland.shp", package="spData")[1], quiet=TRUE)
+auckland <- st_read(system.file("shapes/auckland.gpkg", package="spData")[1], quiet=TRUE)
 res <- probmap(auckland$M77_85, 9*auckland$Und5_81)
 rt <- sum(auckland$M77_85)/sum(9*auckland$Und5_81)
 ppois_pmap <- numeric(length(auckland$Und5_81))
diff --git a/man/skater.Rd b/man/skater.Rd
index 7dbd1237..3d8eb2d7 100644
--- a/man/skater.Rd
+++ b/man/skater.Rd
@@ -72,7 +72,7 @@ skater(edges, data, ncuts, crit, vec.crit, method = c("euclidean",
 \seealso{See Also as \code{\link{mstree}}}
 \examples{
 ### loading data
-bh <- st_read(system.file("etc/shapes/bhicv.shp",
+bh <- st_read(system.file("etc/shapes/bhicv.gpkg",
       package="spdep")[1], quiet=TRUE)
 st_crs(bh) <- "OGC:CRS84"
 ### data standardized 
diff --git a/man/sp.correlogram.Rd b/man/sp.correlogram.Rd
index af15718e..38ba3ff7 100644
--- a/man/sp.correlogram.Rd
+++ b/man/sp.correlogram.Rd
@@ -52,7 +52,7 @@ of the Institute of British Geographers}, 66, 95--118.}
 \seealso{\code{\link{nblag}}, \code{\link{moran}}, \code{\link[stats]{p.adjust}}}
 
 \examples{
-nc.sids <- st_read(system.file("shapes/sids.shp", package="spData")[1], quiet=TRUE)
+nc.sids <- st_read(system.file("shapes/sids.gpkg", package="spData")[1], quiet=TRUE)
 rn <- as.character(nc.sids$FIPS)
 ncCC89_nb <- read.gal(system.file("weights/ncCC89.gal", package="spData")[1],
  region.id=rn)
diff --git a/man/subset.listw.Rd b/man/subset.listw.Rd
index 32c020ce..ae1b573c 100644
--- a/man/subset.listw.Rd
+++ b/man/subset.listw.Rd
@@ -28,7 +28,7 @@ pre <- nb2listw(col.gal.nb)
 print(pre)
 post <- subset(pre, !(1:length(col.gal.nb) \%in\% to.be.dropped))
 print(post)
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 nb <- poly2nb(columbus)
 lw <- nb2listw(nb, style="W")
 attr(lw, "zero.policy")
diff --git a/man/subset.nb.Rd b/man/subset.nb.Rd
index d7379a27..15718be7 100644
--- a/man/subset.nb.Rd
+++ b/man/subset.nb.Rd
@@ -23,7 +23,7 @@
 \seealso{\code{\link{nb2listw}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 coords <- st_coordinates(st_centroid(columbus))
 plot(col.gal.nb, coords)
diff --git a/man/summary.nb.Rd b/man/summary.nb.Rd
index 7bb6fa3b..a4877fac 100644
--- a/man/summary.nb.Rd
+++ b/man/summary.nb.Rd
@@ -35,7 +35,7 @@ the link lengths are also printed. Print and summary functions are also availabl
 \seealso{\code{\link{plot.nb}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 coords <- st_centroid(st_geometry(columbus), of_largest_polygon=TRUE)
 col.gal.nb
diff --git a/man/testnb.Rd b/man/testnb.Rd
index ff665820..6a286abe 100644
--- a/man/testnb.Rd
+++ b/man/testnb.Rd
@@ -41,7 +41,7 @@ A new version of \code{make.sym.nb} by Bjarke Christensen is now included. The o
 \seealso{\code{\link{read.gal}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 col.gal.nb <- read.gal(system.file("weights/columbus.gal", package="spData")[1])
 coords <- st_coordinates(st_centroid(columbus))
 ind <- row.names(as(columbus, "Spatial"))
diff --git a/man/tri2nb.Rd b/man/tri2nb.Rd
index 2f7db98b..76889acc 100644
--- a/man/tri2nb.Rd
+++ b/man/tri2nb.Rd
@@ -38,7 +38,7 @@ The function returns an object of class \code{nb} with a list of integer vectors
 \seealso{\code{\link{knn2nb}}, \code{\link{dnearneigh}}, \code{\link{cell2nb}}}
 
 \examples{
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1], quiet=TRUE)
 coords <- st_centroid(st_geometry(columbus), of_largest_polygon=TRUE)
 ind <- row.names(columbus)
 suppressPackageStartupMessages(require(deldir))
diff --git a/vignettes/CO69.Rmd b/vignettes/CO69.Rmd
index 4bdbac40..1b98ae75 100644
--- a/vignettes/CO69.Rmd
+++ b/vignettes/CO69.Rmd
@@ -50,9 +50,9 @@ run <- require("sp", quiet=TRUE)
 
 ```{r echo=TRUE,eval=run,results='hide'}
 library(spdep)
-eire <- as(sf::st_read(system.file("shapes/eire.shp", package="spData")[1]), "Spatial")
+eire <- as(sf::st_read(system.file("shapes/eire.gpkg", package="spData")[1]), "Spatial")
 row.names(eire) <- as.character(eire$names)
-proj4string(eire) <- CRS("+proj=utm +zone=30 +ellps=airy +units=km")
+#proj4string(eire) <- CRS("+proj=utm +zone=30 +ellps=airy +units=km")
 ```
 ```{r echo=TRUE,eval=run}
 class(eire)
@@ -721,4 +721,4 @@ data sets of this size and shape.
 ## References
 
 [^1]: cropped scans of tables are available from
-    <https://github.com/rsbivand/CO69>.
\ No newline at end of file
+    <https://github.com/rsbivand/CO69>.
diff --git a/vignettes/nb.Rmd b/vignettes/nb.Rmd
index 5eab804d..396aa465 100644
--- a/vignettes/nb.Rmd
+++ b/vignettes/nb.Rmd
@@ -45,7 +45,14 @@ read into R.
 
 ```{r}
 library(spdep)
-NY8 <- as(sf::st_read(system.file("shapes/NY8_utm18.shp", package="spData")), "Spatial")
+if (packageVersion("spData") >= "2.3.2") {
+    NY8a <- sf::st_read(system.file("shapes/NY8_utm18.gpkg", package="spData"))
+} else {
+    NY8a <- sf::st_read(system.file("shapes/NY8_bna_utm18.gpkg", package="spData"))
+    sf::st_crs(NY8a) <- "EPSG:32618"
+    NY8a$Cases <- NY8a$TRACTCAS
+}
+NY8 <- as(NY8a, "Spatial")
 NY_nb <- read.gal(system.file("weights/NY_nb.gal", package="spData"), region.id=as.character(as.integer(row.names(NY8))-1L))
 ```
 
@@ -57,8 +64,8 @@ to the full data set. We retrieve the part of the neighbour list
 in Syracuse using the `subset` method.
 
 ```{r}
-Syracuse <- NY8[NY8$AREANAME == "Syracuse city",]
-Sy0_nb <- subset(NY_nb, NY8$AREANAME == "Syracuse city")
+Syracuse <- NY8[!is.na(NY8$AREANAME) & NY8$AREANAME == "Syracuse city",]
+Sy0_nb <- subset(NY_nb, !is.na(NY8$AREANAME) & NY8$AREANAME == "Syracuse city")
 summary(Sy0_nb)
 ```
 
diff --git a/vignettes/nb_sf.Rmd b/vignettes/nb_sf.Rmd
index de13feb8..ac0cda17 100644
--- a/vignettes/nb_sf.Rmd
+++ b/vignettes/nb_sf.Rmd
@@ -111,7 +111,13 @@ spdep::poly2nb uses two heuristics, first to find candidate neighbours from inte
 Next, we explore a further possible source of differences in neighbour object reproduction, using the original version of the tract boundaries used in ASDAR, but with some invalid geometries as mentioned earlier:
 
 ```{r, echo=dothis, eval=dothis}
-NY8_sf_old <- st_read(system.file("shapes/NY8_utm18.shp", package="spData"), quiet=TRUE)
+if (packageVersion("spData") >= "2.3.2") {
+    NY8_sf_old <- sf::st_read(system.file("shapes/NY8_utm18.gpkg", package="spData"))
+} else {
+    NY8_sf_old <- sf::st_read(system.file("shapes/NY8_bna_utm18.gpkg", package="spData"))
+    sf::st_crs(NY8_sf_old) <- "EPSG:32618"
+    NY8_sf_old$Cases <- NY8_sf_old$TRACTCAS
+}
 table(st_is_valid(NY8_sf_old))
 ```
 
diff --git a/vignettes/sids.Rmd b/vignettes/sids.Rmd
index 6c1eec13..a2777765 100644
--- a/vignettes/sids.Rmd
+++ b/vignettes/sids.Rmd
@@ -44,20 +44,20 @@ library(spdep)
 
 We will be using the **spdep** and **spreg** packages, here version: `r spdep()[1]`, the **sf** package and the **tmap** package. The data from the sources
 referred to above is documented in the help page for the `nc.sids`
-data set in **spData**. The actual data, included in a shapefile of the county boundaries for North Carolina were made available in the **maptools** package [^1]. These data are known to be geographical coordinates (longitude-latitude in decimal degrees) and are assumed to use the NAD27 datum.
+data set in **spData**. The actual data, included in a shapefile of the county boundaries for North Carolina were made available in the **maptools** package [^1]. These data are known to be geographical coordinates (longitude-latitude in decimal degrees) and are assumed to use the NAD27 datum. The ESRI Shapefile is deprecated, and was replaced here by a GeoPackage, written from reading the original files in **spData** 2.3.1:
 
 ```{r echo=TRUE,eval=TRUE}
 library(spdep)
-nc <- st_read(system.file("shapes/sids.shp", package="spData")[1], quiet=TRUE)
-st_crs(nc) <- "EPSG:4267"
+nc <- st_read(system.file("shapes/sids.gpkg", package="spData")[1], quiet=TRUE)
+#st_crs(nc) <- "EPSG:4267"
 row.names(nc) <- as.character(nc$FIPSNO)
 ```
 
-The shapefile format presupposes that you have three files with
+The shapefile format presupposed that you had three files with
 extensions `.shp`, `.shx`, and `.dbf`, where the first
 contains the geometry data, the second the spatial index, and the third
-the attribute data. They are required to have the same name apart from
-the extension, and are read here using `sf::st_read()` into the
+the attribute data. They were required to have the same name apart from
+the extension, and were read here using `sf::st_read()` into the
 `sf` object `nc`; the class is defined in **sf**. The centroids of the largest polygon in each county are available
 using the `st_centroid` method from **sf** as an **sfc** POINT object, and can
 be used to place labels after the extraction of the coordinate matrix:
@@ -463,4 +463,4 @@ tables of covariates rather than geographical grids.
 ## References
 
 [^1]: These data were taken with permission from a now-offline link:
-[sal.agecon.uiuc.edu/datasets/sids.zip]; see also [GeoDa Center](https://geodacenter.github.io/data-and-lab/) for a contemporary source.
\ No newline at end of file
+[sal.agecon.uiuc.edu/datasets/sids.zip]; see also [GeoDa Center](https://geodacenter.github.io/data-and-lab/) for a contemporary source.