mstrimas / smoothr

Spatial Feature Smoothing in R
http://strimas.com/smoothr
GNU General Public License v3.0
96 stars 5 forks source link

Invalid geometry from ksmooth #1

Open andrew-plowright opened 5 years ago

andrew-plowright commented 5 years ago

Hello,

I tried using the ksmooth smoothing function on a SHP file (uploaded to Google Drive, follow link below). I checked to ensure that the input (x) is valid. However, the output (y) contains invalid geometry, and cannot be saved back into a SHP file.

https://drive.google.com/file/d/1w7P-oRr7EohAcX0nE71zQza53nmjUb7f/view?usp=sharing

Thanks!

library(sf)
library(smoothr)

x <- st_read("makes_invalid2.shp")

any(!st_is_valid(x))

y <- smooth(x, method = "ksmooth", smoothness = 4)

any(!st_is_valid(y))
mstrimas commented 5 years ago

Hi Andrew, thanks for the bug report! I'm not too surprised since I'm using a fairly simple approach for the smoothing and I don't have a great grasp about the how these sf objects are stored or about geometric operations on then.

Your original shapefile has several geometries each composed of multiple polygons. As a starting point, I was able to extract just the one polygon that's causing the problem.

library(sf)
library(smoothr)

x <- st_read("~/Downloads/makes_invalid2/makes_invalid2.shp")[21,] %>% 
  st_geometry()
x <- st_polygon(list(x[[1]][[1]]))

any(!st_is_valid(x))

y <- smooth(x, method = "ksmooth", smoothness = 4)

any(!st_is_valid(y))

We can also get a reason for the invalidity:

st_is_valid(y, reason = TRUE)
#> "Invalid Coordinate[nan nan]"

But that's as far as I've gotten, I still can't pinpoint exactly what the issue is. I'll keep fiddling around, although it's possible this is over my head.

andrew-plowright commented 5 years ago

Thanks for looking into it!

It seems like a sequence of 428 vertices are invalidated by the smoothing process. Not sure why, though.

x <- st_read("C:/Users/Andrew/Desktop/scratch/makes_invalid2.shp")[21,] %>% 
  st_geometry()
x <- st_polygon(list(x[[1]][[1]]))

y <- smooth(x, method = "ksmooth", smoothness = 4)

# Determine which vertices are invalid
invalid <- apply(y[[1]],1 , function(x) any(is.na(x)))