diff --git a/DESCRIPTION b/DESCRIPTION
index d1ab828..6207cf9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: spatialreg
-Version: 1.3-3
-Date: 2024-05-31
+Version: 1.3-4
+Date: 2024-06-10
 Title: Spatial Regression Analysis
 Encoding: UTF-8
 Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"), email = "Roger.Bivand@nhh.no", comment=c(ORCID="0000-0003-2392-6140")),
@@ -25,7 +25,7 @@ Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"), email = "Roger.Bi
              person(given = "Mauricio", family = "Sarrias", role = c("ctb"), email = "mauricio.sarrias@ucn.cl"),
              person(given = "JuanTomas", family = "Sayago", role = c("ctb"), email = "juantomas.sayago@gmail.com"),
              person("Michael", "Tiefelsdorf", role = "ctb"))
-Depends: R (>= 3.3.0), spData, Matrix, sf
+Depends: R (>= 3.3.0), spData (>= 2.3.1), Matrix, sf
 Imports: spdep (>= 1.3-1), coda, methods, MASS, boot, splines, LearnBayes,
         nlme, multcomp
 Suggests: parallel, RSpectra, tmap, foreign, spam, knitr, lmtest, expm,
diff --git a/NEWS.md b/NEWS.md
index 0214ed7..5510430 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,8 @@
-# Version 1.3-3 (development)
+# Version 1.3-4 (development)
+
+* migrate from ESRI Shapefile to GeoPackage #50
+
+# Version 1.3-3 (2024-05-31)
 
 * protect `errorsarlm` against missing `Durbin=` if only intercept
 
diff --git a/man/ME.Rd b/man/ME.Rd
index 5de00f8..5d5eb4d 100644
--- a/man/ME.Rd
+++ b/man/ME.Rd
@@ -65,7 +65,7 @@ coef(summary(glmME))
 anova(glmbase, glmME, test="Chisq")
 \dontrun{
 require("sf", quietly=TRUE)
-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)
 #require("spdep", quietly=TRUE)
 col.gal.nb <- spdep::read.gal(system.file("weights/columbus.gal", package="spData")[1])
 lw <- spdep::nb2listw(col.gal.nb)
@@ -96,7 +96,7 @@ COL.ME.NA <- ME(CRIME ~ INC + HOVAL, data=NA.columbus, family="gaussian",
 COL.ME.NA$na.action
 summary(lm(CRIME ~ INC + HOVAL + fitted(COL.ME.NA), data=NA.columbus,
  na.action=na.exclude))
-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 <- spdep::read.gal(system.file("weights/ncCC89.gal", package="spData")[1],
  region.id=rn)
diff --git a/man/SpatialFiltering.Rd b/man/SpatialFiltering.Rd
index 0081984..f1aa340 100644
--- a/man/SpatialFiltering.Rd
+++ b/man/SpatialFiltering.Rd
@@ -55,7 +55,7 @@ SpatialFiltering(formula, lagformula=NULL, data=list(), na.action=na.fail,
 \seealso{\code{\link[stats]{lm}}, \code{\link[base]{eigen}}, \code{\link{nb2listw}}, \code{\link{listw2U}}}
 \examples{
 require("sf", quietly=TRUE)
-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)
 #require("spdep", quietly=TRUE)
 col.gal.nb <- spdep::read.gal(system.file("weights/columbus.gal", package="spData")[1])
 lmbase <- lm(CRIME ~ INC + HOVAL, data=columbus)
diff --git a/man/aple.Rd b/man/aple.Rd
index 314b18c..f32ebd5 100644
--- a/man/aple.Rd
+++ b/man/aple.Rd
@@ -27,7 +27,7 @@ aple(x, listw, override_similarity_check=FALSE, useTrace=TRUE)
 \seealso{\code{\link[spdep]{nb2listw}}, \code{\link{aple.mc}}, \code{\link{aple.plot}}
 }
 \examples{
-wheat <- st_read(system.file("shapes/wheat.shp", package="spData")[1], quiet=TRUE)
+wheat <- st_read(system.file("shapes/wheat.gpkg", package="spData")[1], quiet=TRUE)
 library(spdep)
 nbr1 <- spdep::poly2nb(wheat, queen=FALSE)
 nbrl <- spdep::nblag(nbr1, 2)
diff --git a/man/aple.mc.Rd b/man/aple.mc.Rd
index cb533c1..ec9fab4 100644
--- a/man/aple.mc.Rd
+++ b/man/aple.mc.Rd
@@ -26,7 +26,7 @@ aple.mc(x, listw, nsim, override_similarity_check=FALSE, useTrace=TRUE)
 \seealso{\code{\link{aple}}, \code{\link[boot]{boot}}}
 \examples{
 \dontrun{
-wheat <- st_read(system.file("shapes/wheat.shp", package="spData")[1], quiet=TRUE)
+wheat <- st_read(system.file("shapes/wheat.gpkg", package="spData")[1], quiet=TRUE)
 nbr1 <- spdep::poly2nb(wheat, queen=FALSE)
 nbrl <- spdep::nblag(nbr1, 2)
 nbr12 <- spdep::nblag_cumul(nbrl)
diff --git a/man/aple.plot.Rd b/man/aple.plot.Rd
index 7f96bfb..70914aa 100644
--- a/man/aple.plot.Rd
+++ b/man/aple.plot.Rd
@@ -35,7 +35,7 @@ localAple(x, listw, override_similarity_check=FALSE, useTrace=TRUE)
 \seealso{\code{\link{aple}}}
 \examples{
 \dontrun{
-wheat <- st_read(system.file("shapes/wheat.shp", package="spData")[1], quiet=TRUE)
+wheat <- st_read(system.file("shapes/wheat.gpkg", package="spData")[1], quiet=TRUE)
 nbr1 <- spdep::poly2nb(wheat, queen=FALSE)
 nbrl <- spdep::nblag(nbr1, 2)
 nbr12 <- spdep::nblag_cumul(nbrl)
diff --git a/man/impacts.Rd b/man/impacts.Rd
index 858275b..06628d1 100644
--- a/man/impacts.Rd
+++ b/man/impacts.Rd
@@ -95,7 +95,7 @@ Roger Bivand, Gianfranco Piras (2015). Comparing Implementations of Estimation M
 }
 \examples{
 require("sf", quietly=TRUE)
-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)
 #require("spdep", quietly=TRUE)
 col.gal.nb <- spdep::read.gal(system.file("weights/columbus.gal", package="spData")[1])
 listw <- spdep::nb2listw(col.gal.nb)
diff --git a/man/sarlm_tests.Rd b/man/sarlm_tests.Rd
index 2681193..e9d053b 100644
--- a/man/sarlm_tests.Rd
+++ b/man/sarlm_tests.Rd
@@ -76,7 +76,7 @@ Dordrecht: Kluwer, pp. 121--122.
 
 \examples{
 require("sf", quietly=TRUE)
-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)
 #require("spdep", quietly=TRUE)
 col.gal.nb <- spdep::read.gal(system.file("weights/columbus.gal", package="spData")[1])
 lm.mod <- lm(CRIME ~ HOVAL + INC, data=columbus)
diff --git a/man/sparse_mat.Rd b/man/sparse_mat.Rd
index 1e41927..738a511 100644
--- a/man/sparse_mat.Rd
+++ b/man/sparse_mat.Rd
@@ -46,7 +46,7 @@ powerWeights(W, rho, order=250, X, tol=.Machine$double.eps^(3/5))
 \examples{
 \dontrun{
 require(sf, quietly=TRUE)
-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)
 #require(spdep, quietly=TRUE)
 col.gal.nb <- spdep::read.gal(system.file("weights/columbus.gal", package="spData")[1])
 col.listw <- spdep::nb2listw(col.gal.nb)
diff --git a/man/spautolm.Rd b/man/spautolm.Rd
index 60f4b4f..0d675a6 100644
--- a/man/spautolm.Rd
+++ b/man/spautolm.Rd
@@ -199,7 +199,7 @@ summary(ecar1wM)
 }
 \dontrun{
 require("sf", quietly=TRUE)
-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)
 ft.SID74 <- sqrt(1000)*(sqrt(nc.sids$SID74/nc.sids$BIR74) +
  sqrt((nc.sids$SID74+1)/nc.sids$BIR74))
 lm_nc <- lm(ft.SID74 ~ 1)
diff --git a/man/trW.Rd b/man/trW.Rd
index 0ce5798..5243f32 100644
--- a/man/trW.Rd
+++ b/man/trW.Rd
@@ -38,7 +38,7 @@ A numeric vector of \code{m} traces, with \dQuote{timings} and \dQuote{type} att
 }
 \examples{
 require("sf", quietly=TRUE) 
-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)
 #require(spdep, quietly=TRUE)
 col.gal.nb <- spdep::read.gal(system.file("weights/columbus.gal", package="spData")[1])
 listw <- spdep::nb2listw(col.gal.nb)
diff --git a/vignettes/SpatialFiltering.Rmd b/vignettes/SpatialFiltering.Rmd
index 15c0ca9..11d9534 100644
--- a/vignettes/SpatialFiltering.Rmd
+++ b/vignettes/SpatialFiltering.Rmd
@@ -54,7 +54,13 @@ by Pedro Peres-Neto and is discussed in @dray+legendre+peres-neto:06 and
 ```{r echo=TRUE}
 library(spdep)
 require("sf", quietly=TRUE)
-NY8 <- st_read(system.file("shapes/NY8_utm18.shp", package="spData"))
+if (packageVersion("spData") >= "2.3.2") {
+    NY8 <- sf::st_read(system.file("shapes/NY8_utm18.gpkg", package="spData"))
+} else {
+    NY8 <- sf::st_read(system.file("shapes/NY8_bna_utm18.gpkg", package="spData"))
+    sf::st_crs(NY8) <- "EPSG:32618"
+    NY8$Cases <- NY8$TRACTCAS
+}
 NY_nb <- read.gal(system.file("weights/NY_nb.gal", package="spData"), override.id=TRUE)
 ```
 
@@ -89,22 +95,18 @@ set.seed(111)
 ```
 
 ```{r, echo=TRUE}
-nyME <- ME(Cases~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8, offset=log(POP8), family="poisson", listw=NYlistwW, alpha=0.44)
+nyME <- ME(Cases~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8, offset=log(POP8), family="poisson", listw=NYlistwW, alpha=0.46)
 ```
 
 ```{r, echo=TRUE}
 nyME
-NY8$eigen_24 <- fitted(nyME)[,1]
-NY8$eigen_223 <- fitted(nyME)[,2]
+NY8$eigen_1 <- fitted(nyME)[,1]
+NY8$eigen_2 <- fitted(nyME)[,2]
 ```
 
-```{r, results='asis',echo=FALSE, eval=FALSE}
-.iwidth <- 6
-.iheight <- 4
-.ipointsize <- 10
-library(RColorBrewer)
+```{r, echo=TRUE}
 #gry <- brewer.pal(9, "Greys")[-1]
-spplot(NY8, c("eigen_24", "eigen_223"), col.regions=grey.colors(6, 0.95, 0.55, 2.2), cuts=5)
+plot(NY8[,c("eigen_1", "eigen_2")])
 ``` 
 
 ```{r, echo=TRUE}
diff --git a/vignettes/nb_igraph.Rmd b/vignettes/nb_igraph.Rmd
index bd680a5..c1dbcf8 100644
--- a/vignettes/nb_igraph.Rmd
+++ b/vignettes/nb_igraph.Rmd
@@ -50,7 +50,7 @@ if (dothis) {
 ```
 ```{r, echo=dothis, eval=dothis}
 library(sf)
-columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1])
+columbus <- st_read(system.file("shapes/columbus.gpkg", package="spData")[1])
 row.names(columbus)[1:10]
 ```
 
diff --git a/vignettes/sids_models.Rmd b/vignettes/sids_models.Rmd
index a431c80..d975482 100644
--- a/vignettes/sids_models.Rmd
+++ b/vignettes/sids_models.Rmd
@@ -49,8 +49,8 @@ data set in **spData**. The actual data, included in a shapefile of the county b
 
 ```{r echo=TRUE,eval=TRUE}
 library(spdep)
-nc <- st_read(system.file("shapes/sids.shp", package="spData")[1], quiet=TRUE)
-st_crs(nc) <- "+proj=longlat +datum=NAD27"
+nc <- st_read(system.file("shapes/sids.gpkg", package="spData")[1], quiet=TRUE)
+#st_crs(nc) <- "+proj=longlat +datum=NAD27"
 row.names(nc) <- as.character(nc$FIPSNO)
 ```
 
@@ -175,4 +175,4 @@ plot(ll ~ lambda, ecarIawll$llprof, type="l")
 ## 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.