Message-ID: <CAL_d4CW5_Mgk_16FWppAcP8UExt6TDesPxY1+4aigL6mj0gtNQ@mail.gmail.com>
Date: 2025-12-19T06:15:10Z
From: Aymard Loïc KWEKEU KWEKEU
Subject: Warnings and error during the creation of the neighborhood matrix by recoding old code using deprecated rgdal in CRAN with sf and sp
In-Reply-To: <e14b278a-1bfc-7b4b-a531-77cf62792776@reclus2.nhh.no>
Hello everyone, sorry for my late reply, my supervisor and I have
determined that this was just an isolated observation and we have decided
to delete it.
Le jeu. 18 d?c. 2025, 09:52, Roger Bivand <Roger.Bivand at nhh.no> a ?crit :
> On Thu, 18 Dec 2025, Edzer Pebesma wrote:
>
> > I think the CRAN package sfislands tries to help solve the issue of
> > unconnected graphs:
>
> It does - see also https://journal.r-project.org/articles/RJ-2025-015/,
> so
> does the spdep vignette using distances between polygon boundaries, but
> the primary cause here may be the differences between rgdal/sp with an
> unknown version of GDAL (most likely much older) and sf also with an
> unknown version of GDAL (almost certainly later) in reading the
> coordinates from the ESRI shapefile.
>
> Without access to the shapefile, I don't think that we can conclude that
> the resolution is not rather to use sf with a properly defined CRS (snap
> defaults differ a little in spdep now for planar and spherical
> coordinates), and examine the impact of slackening snap by a couple of
> centimetres.
>
> Roger
>
> >
> > https://cran.r-project.org/web/packages/sfislands/index.html
> >
> > On 17/12/2025 19:42, Josiah Parry wrote:
> >> 200% echo everything said by Dr. Bivand below.
> >>
> >> In addition, I would focus on the error message:
> >>
> >> Empty neighbor sets found (zero.policy: FALSE)
> >>
> >>
> >> This is the thing you need to solve. Why are there locations without
> >> neighbors? A "cheap way out" which we offer in ArcGIS Pro is to
> provide a
> >> union of the neighbor sets using KNN and contiguity. This ensures that
> >> each
> >> location has at minimum k number of neighbors. You can accomplish this
> >> like
> >> the example below. *However*, I would consider the implications of
> using
> >> distance metrics for polygons. See if what you're doing actually makes
> >> sense conceptually.
> >>
> >> library(spdep)
> >>
> >> geo <- system.file("shapes/columbus.gpkg", package="spData") |>
> >> sf:: st_read(quiet = TRUE) |>
> >> sf:: st_geometry()
> >>
> >> nb_contig <- poly2nb(geo)
> >>
> >> # helper from sfdep for st_centroid() -> knearneigh() -> knn2nb()
> >> nb_k1 <- sfdep::st_knn(geo)
> >>
> >> union.nb(nb_contig, nb_k1)
> >>
> >> On Wed, Dec 17, 2025 at 1:06?AM Roger Bivand via R-sig-Geo <
> >> r-sig-geo at r-project.org> wrote:
> >>
> >>> On Tue, 16 Dec 2025, Aymard Lo?c KWEKEU KWEKEU wrote:
> >>>
> >>> First, you do not show the versions of software being used, including
> the
> >>> platform, versions of rgdal, sp, sf, spdep, mgcv, R2BayesX, possibly
> >>> others. The versions of GDAL used in rgdal and sf may also matter, as
> >>> GDAL
> >>> also changes over time (the precision of coordinates in the same ESRI
> >>> Shapefile may differ, artefacts in polygons bay be corrected
> >>> differently).
> >>>
> >>> Comments inline.
> >>>
> >>>>
> >>>> *Initial code with rgdal's readOGR function: *
> >>>>
> >>>> geo <- readOGR("./Data/Carto/COMMUNE_CARTO.shp")
> >>>
> >>> Without an example file (all the files making up the ESRI shapefile,
> at
> >>> least *.shp, *.shx, *.dbf, preferably with *.prj), there is no way to
> >>> establish anything.
> >>>
> >>>>
> >>>> geo2 <- subset(geo, geo at data$INSEE_COM %in% selected_com)
> >>>>
> >>>
> >>> Never refer to S4 slots using @ in workflows - geo$INSEE_COM has been
> >>> sufficient now for 20 years.
> >>>
> >>>> DATA2BND <- sp2bnd(geo2, geo2 at data$INSEE_COM)
> >>>>
> >>>> plot(DATA2BND)
> >>>>
> >>>> list_names <- as.data.frame(names(DATA2BND)) %>%
> >>>> group_by(names(DATA2BND)) %>%
> >>>> summarise(n = n())
> >>>>
> >>>
> >>> sp2bnd is a function in R2BayesX and is strictly superfluous to this
> >>> question.
> >>>
> >>>> nb <- mgcv:::pol2nb(DATA2BND)$nb
> >>>>
> >>>
> >>> mgcv:::pol2nb is a function in mgcv and is strictly superfluous to
> this
> >>> question.
> >>>
> >>>> W.nb <- poly2nb(
> >>>> geo2,
> >>>> row.names = rownames(geo2 at data$INSEE_COM)
> >>>> )
> >>>>
> >>>> w <- nb2mat(W.nb, style = "B")
> >>>
> >>> This is the core of the question - why no warnings or errors?
> >>>
> >>>>
> >>>>
> >>>> *New code with sf and sp: *
> >>>>
> >>>> geo_sf <- st_read("./Data/Carto/COMMUNE_CARTO.shp", quiet = TRUE)
> >>>
> >>> Coercing to an sp-class is superfluous.
> >>>
> >>>> geo <- as(geo_sf, "Spatial")
> >>>> geo2 <- subset(geo, geo at data$INSEE_COM %in% selected_com)
> >>>>
> >>>> DATA2BND <- sp2bnd(geo2, geo2 at data$INSEE_COM)
> >>>>
> >>>> plot(DATA2BND)
> >>>>
> >>>> list_names <- as.data.frame(names(DATA2BND)) %>%
> >>>> group_by(names(DATA2BND)) %>%
> >>>> summarise(n = n())
> >>>>
> >>>> nb <- mgcv:::pol2nb(DATA2BND)$nb
> >>>
> >>> It is possible that the coordinate values of the polygons in geo2
> differ
> >>> here by a very small amount if the versions of GDAL used to install
> rgdal
> >>> and sf differ, and also because sf's implementation for reading
> >>> geometries
> >>> is more modern and cleaner than that in rgdal (>20 years old).
> >>>
> >>>>
> >>>> W.nb <- poly2nb(geo2, row.names = rownames(geo2 at data$INSEE_COM),
> >>> queen=T)
> >>>> Warning messages:
> >>>>
> >>>> 1: In poly2nb(geo2, row.names = rownames(geo2 at data$INSEE_COM),
> queen =
> >>> T):
> >>>>
> >>>> some observations have no neighbors;
> >>>>
> >>>> if this seems unexpected, try increasing the snap argument.
> >>>
> >>> Recent changes in spdep warn about observations with no neighbours,
> and
> >>> about subgraphs. Please read the relevant vignette:
> >>>
> https://cran.r-project.org/web/packages/spdep/vignettes/subgraphs.html
> >>>
> >>> No-neighbour observations may be detected when apparently contiguous
> >>> polygons are separated by >10mm (default value) if the input polygons
> are
> >>> an sf object, otherwise sqrt(.Machine$double.eps), (so the sp case is
> >>> unchanged).
> >>>
> >>>>
> >>>> 2: In poly2nb(geo2, row.names = rownames(geo2 at data$INSEE_COM),
> queen =
> >>> T):
> >>>>
> >>>> neighbor object has 7 sub-graphs;
> >>>>
> >>>> if this sub-graph count seems unexpected, try increasing the snap
> >>> argument.
> >>>>
> >>>> w_SMD_2025 <- nb2mat(W.nb, style = "B")
> >>>>
> >>>> Error in nb2listw(neighbors, glist = glist, style = style,
> zero.policy =
> >>>> zero.policy):
> >>>>
> >>>> Empty neighbor sets found (zero.policy: FALSE)
> >>>
> >>> If you accept a neighbour object with subgraphs (probably no-neighbour
> >>> observations), set zero.policy=TRUE in nb2mat. This may however cause
> >>> problems later on.
> >>>
> >>> I suggest rather avoiding sp, subsetting geo_sf, and using that
> object in
> >>> poly2nb. You also need to check whether the polygons are planar or
> >>> spherical, as this affects the default value of snap.
> >>>
> >>> If you need to share data, do not attach to posts on this list, but
> >>> rather
> >>> raise an issue on
> >>> https://github.com/r-spatial/spdep/issues/,
> >>> where the
> >>> data and code for a minimal reproducible example may be attached.
> >>>
> >>> Hope this clarifies,
> >>>
> >>> Roger
> >>>
> >>> --
> >>> Roger Bivand
> >>> Emeritus Professor
> >>> Department of Economics, Norwegian School of Economics,
> >>> Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
> >>> e-mail: Roger.Bivand at nhh.no
> >>> _______________________________________________
> >>> R-sig-Geo mailing list
> >>> R-sig-Geo at r-project.org
> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>>
> >>
> >> [[alternative HTML version deleted]]
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo at r-project.org
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> >
>
> --
> Roger Bivand
> Emeritus Professor
> Department of Economics, Norwegian School of Economics,
> Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
> e-mail: Roger.Bivand at nhh.no
>
[[alternative HTML version deleted]]