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