Unable to snap points to linestring in R due to sf package error

Question

I have two shapefiles, points and linestring, as shown below. In my R (version 4.2 on Windows 10), I am trying to snap all the points to the nearest line and getting error.

library(maptools)
library(rgdal)
library(sp)
library(sf)

# setwd(path_to_dir)
points <- read_sf("points.shp")
line <- read_sf("line.shp")
# Both files are in the same UTM projection CRS

# Drop ZM
points$geometry <- st_zm(points$geometry)

# Snap
pts_snap <- snapPointsToLines(points, line, maxDist=50, withAttrs = TRUE, idField="IDField")

I get the following error:

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘is.projected’ for signature ‘"sf"’

I found a relevant question here, and tried steps suggested there to fix the error, but no luck.

class(points)
# [1] "sf"         "tbl_df"     "tbl"        "data.frame"
class(line)
# [1] "sf"         "tbl_df"     "tbl"        "data.frame"

# Transform to WGS 84
points <- st_transform(points,4326)
line <-   st_transform(line,4326)

shape_proj1 <-st_transform(points, CRS("+proj=gnom +lat_0=90 +lon_0=-50"))
shape_proj2 <-st_transform(line, CRS("+proj=gnom +lat_0=90 +lon_0=-50"))

# and tried Snapping again, but the same error

Can anyone help me to fix the error?

Sample scenario:

Edit 1:

Here is ST_CRS(layer) info.

> sf::st_crs(points)
Coordinate Reference System:
  User input: ETRS89 / UTM zone 32N 
  wkt:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        DATUM["European Terrestrial Reference System 1989",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[......

What were `sf::st_crs(point)` and `sf::st_crs(line)` immediately after `read_sf(`? — Chris, Oct 29 '22 at 16:48
@Chris: Just updated the question with `st_crs(points)` info. It's the same for line. — timetraveller007, Oct 29 '22 at 17:04
Any chance you can share the shapefiles? Or a least `dput(head(points))`, `dput(line)`? — Grzegorz Sapijaszko, Oct 29 '22 at 17:15
can you try it in a new R session, but without loading `sp`, `rgdal` or `maptools`, and use `sf::st_snap()`? — tospig, Oct 30 '22 at 10:15
@tospig: +1, mate! It worked, but didn't solve it completely. When I use `sf` package in a separate R session and `st_snap(points, line, 30)`. However, the snapped points are only at the start and the end point of the line. No points in between, thus I am not sure if all the points are snapped to start/end of the line... — timetraveller007, Oct 30 '22 at 11:05
Ok, try now by replacing `sf::st_snap()` with `maptools::snapPointsToLine` (where you don't call `library(maptools)`, but rather explicitly use `maptools::` to call the function — tospig, Oct 31 '22 at 00:44

dimfalk · Accepted Answer · 2022-10-31T02:59:59.067

Let's make this reproducible first:

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE

set.seed(42)

mask <- data.frame(lon = c(8, 15),
                   lat = c(52, 63)) |> 
  st_as_sf(coords = c("lon", "lat"), crs = "epsg:4326") |> 
  st_bbox() |> 
  st_as_sfc() |> 
  st_transform("epsg:25832")

poly <- data.frame(lon = c(10, 12),
                   lat = c(55, 60)) |> 
  st_as_sf(coords = c("lon", "lat"), crs = "epsg:4326") |> 
  st_bbox() |> 
  st_as_sfc() |> 
  st_transform("epsg:25832")

n <- 100

points <- st_sample(poly, n)

line <- data.frame(lon = c(13, 13),
                   lat = c(54, 61)) |> 
  st_as_sf(coords = c("lon", "lat"), crs = "epsg:4326") |> 
  st_combine() |> 
  st_cast("LINESTRING") |> 
  st_transform("epsg:25832")

plot(mask)
plot(points, add = TRUE)
plot(line, col = "red", add = TRUE)

As already mentioned, st_snaps() only seems to snap to the start/end nodes of the linestring. Still not sure if this correct, but a workaround for this exists here.

# determine nearest points and cast the lines to points
nrst <- st_nearest_points(points, line) |> 
  st_cast("POINT")

# keep only the end nodes of the connecting lines
p_snapped <- nrst[seq(2, n * 2, 2)]

plot(p_snapped, col = "blue", add = TRUE)

Unable to snap points to linestring in R due to sf package error

1 Answers1