Search code examples
rvectorgeometryterra

Problems working with vector layers in R terra "Points of LinearRing do not form a closed linestring "

in the last weeks I got inconsistent results using vector layers in R using the terra package. For example, using this dataset: http://www.soest.hawaii.edu/pwessel/gshhg/gshhg-shp-2.3.7.zip

A month ago I performed these actions (I used GSHHS_f_L1.shp but this problem happened with several layers in the zip file):

library(terra)
xx<-vect('~/GSHHS_f_L1.shp')
plot(xx)
is.valid(xx)
TRUE

The plot and everything else worked fine enter image description here

One week ago I did it again using the same layer. The plot worked fine, but

is.valid(xx)
FALSE

yy<-makeValid(xx)
is.valid(yy)
TRUE

and again, everything worked fine. However, the last days I got this error (same layer) that I still can't solved

is.valid(xx)
Error: IllegalArgumentException: Points of LinearRing do not form a closed linestring

yy<-makeValid(xx)
Error: IllegalArgumentException: Points of LinearRing do not form a closed linestring

After this, the plot is the only action that works. I try the common solution of creating a 0-buffer, but I got the same error. I updated terra and downloaded the zip file again, but nothing change. I know of some alternatives using sf, but I feel very comfortable with terra and I want to know if there is a solution using this package. I'm using R 4.4.0, terra 1.7.78 in a laptop Linux Ubuntu 22.04.4 LTS. The only change in the last weeks was that I installed the package gower from github. Any recommendation will be very appreciated

Solution

  • On my Win10 machine, I found something weird going on with the data (or perhaps it is a bug?) that occurred when applying an sf solution, but not when using only terra. When I applied terra::makeValid() to the entire dataset, it worked. However, the equivalent sf::st_make_valid() did not work on the full dataset. I have added a potential terra-only workaround inspired by the sf issue that will hopefully work for you on your machine. I have also included the sf option by way of illustration, mainly to highlight the weird behaviour, which may be indicative of the issue you face.

    First, using only terra:

    library(terra)

# Load data
xx <- vect('~/GSHHS_f_L1.shp')

# Check validity
unique(is.valid(xx))
# [1]  TRUE FALSE

# Return valid
xx_is <- xx[is.valid(xx)]

unique(is.valid(xx_is))
# [1] TRUE

# Return invalid
xx_not <- xx[!is.valid(xx)]

is.valid(xx_not)
# [1] FALSE

# Make valid
xx_is1 <- makeValid(xx_not)

is.valid(xx_is1)
# [1] TRUE

# Recreate xx
xx <- vect(c(xx_is, xx_is1))

unique(is.valid(xx))
# [1] TRUE

    Using sf:

    library(terra)
library(sf)
library(dplyr)

xx <- vect('~/GSHHS_f_L1.shp')

# xx SpatVector to sf
xx_sf <- st_as_sf(xx) %>%
  st_make_valid() %>% # Does not work, but subsequent steps do not work without it either
  mutate(valid = st_is_valid(.))

# Check for invalid geometries (should not return any FALSE but does for some reason)
unique(st_is_valid(xx_sf))
# [1]  TRUE FALSE

# Return invalid polygon
not_valid <- filter(xx_sf, !valid)

st_is_valid(not_valid)
# [1] FALSE

# Make valid
is_valid <- st_make_valid(not_valid)

st_is_valid(is_valid)
# [1] TRUE

# Find index of oringinal invalid polygon
intersects <- st_intersects(is_valid, xx_sf)
# Sparse geometry binary predicate list of length 1, where the predicate was
# `intersects'
#  1: 2245

# Return invalid polygon
xx_sf_inv <- xx_sf[2245,]

xx_sf_inv
# Simple feature collection with 1 feature and 7 fields
# Geometry type: MULTIPOLYGON
# Dimension:     XY
# Bounding box:  xmin: -69.78281 ymin: 43.75223 xmax: -69.70734 ymax: 43.88222
# Geodetic CRS:  WGS 84
# id level source parent_id sibling_id     area                       geometry valid
# 2245 2380     1    WVS        -1         -1 45.04373 MULTIPOLYGON (((-69.76185 4... FALSE

# Filter out invalid polygon from xx_sf
xx_sf <- filter(xx_sf, id != 2380)

# Add valid version of polygon back to xx_sf
xx_sf <- rbind(xx_sf, is_valid)

# Recheck geometries
unique(st_is_valid(xx_sf))
# [1] TRUE

# Convert back to SpatVector
xx <- vect(xx_sf)

# Recheck geometries again
unique(is.valid(xx))
# [1] TRUE