Skip to content

Using sf package to calculate line midpoints from Spatial Line data (both LINESTRING and MULTILINESTRING)

3 messages · Julian Parnelle, @bo m@iii@g oii eirick@de, Roger Bivand

Original
Julian Parnelle
# 
Dear list,

My question is straightforward: how could one use the sf package to calculate the Longitude and Latitude of the line midpoints in spatial line data, in a way that works both for LINESTRING and MULTILINESTRING sf data.frames? (one per spatial line, in case the spatial line is composed of multiple line segments). After searching online for quite a bit, the only direct reference I found was this GIS Stack Exchange question from three years ago: https://gis.stackexchange.com/questions/277219/sf-equivalent-of-r-maptools-packages-spatiallinesmidpoints. However, the solution proposed there is only more or less reliable for LINESTRING data and I lack sufficient knowledge of sf objects and MULTILINESTRING representation to update that and make it more reliable.

I was used to doing that with the SpatialLinesMidPoints command, but I am finding it quite difficult to accomplish the same with the sf package. I do know that this task could perhaps be achieved by interpolating across the lines up to 50% of its extension, using for example ST_Line_Interpolate_Point. However, I don't see to be able to achieve exactly the same with either st_line_sample or st_segmentize.

Any suggestions would be very welcome.

Julian
@bo m@iii@g oii eirick@de
# 
Dear Julian,

To get the midpoint of a line data just use st_centroid(). Here is a  
brief example:

library(sf)

# example with two lines
ds <- data.frame(id = 1:2,
              geom = c('LINESTRING(0 0, 2 2)', 'LINESTRING(3 2, 4 5)'))
ds <- st_as_sf(ds, wkt = "geom")
st_centroid(ds)

# example with a multi line string
ds <- data.frame(id = 1:1,
              geom = c('MULTILINESTRING ((0 0, 2 2),(3 2, 4 5))'))
ds <- st_as_sf(ds, wkt = "geom")
st_centroid(ds)

Hope that helps,
Tim

Zitat von Julian Parnelle <julianparnelle at outlook.com>:
Roger Bivand
#
On Wed, 3 Feb 2021, abo at elrick.de wrote:

            

      
      
      
    

        
Not quite, but checking why not and in which cases would be interesting - 
please do run the code, the ouput is what is interesting:

library(maptools)
xx <- readShapeLines(system.file("shapes/fylk-val.shp",
   package="maptools")[1],
   proj4string=CRS("+proj=utm +zone=33 +datum=WGS84"))
spdf <- SpatialLinesMidPoints(xx)
oo <- sf::st_centroid(sf::st_as_sf(xx))
all.equal(sf::st_coordinates(oo), coordinates(spdf),
   check.attributes=FALSE, scale=1)
df <- sf::st_coordinates(oo) - coordinates(spdf)
df[apply(df, 1, function(x) any(abs(x) > 1)),]

geos <- rgeos::gCentroid(xx, byid=TRUE)
all.equal(sf::st_coordinates(oo), coordinates(geos),
   check.attributes=FALSE, scale=1) # same as st_centroid()

And:

lr <- t(sapply(1:length(xx), function(i)
   coordinates(rgeos::gInterpolate(xx[i,], d=0.5, normalized=TRUE))))
all.equal(coordinates(spdf), lr, check.attributes=FALSE, scale=1)

So rgeos::gInterpolate() and SpatialLinesMidPoints() appear to be 
interpolating along the line segments, that is using linear reference.

res <- gDistance(geos, xx, byid=TRUE)
summary(c(res))

res1 <- gDistance(spdf, xx, byid=TRUE)
summary(diag(res1))

So the linear reference measures give points on the lines, but centroids 
do not necessarily do so. Neither of GEOSInterpolateNormalized_r() nor 
GEOSInterpolate_r() appear in code in sf/src, so maptools or rgeos are at 
present the only alternatives; maybe one of the spatstat packages also 
provides linear reference.

Hope this clarifies,

Roger