TransWikia.com

Raster to Lines fails (R)

Geographic Information Systems Asked by Heymans on February 26, 2021

I have the following set of images where a set of 8-neighbor contiguous pixels are outlined in each image, representing a line. The images are distinct, but the lines generally run from left to right. I want to get a single Line (Linestring) that runs through the centers of the pixels, and export as a shape or geoJson file. We imposed an unprojected WGS84 coordinate system (EPSG:4326) on the matrix because we use a geopackage downstream that’ll require it. My question is how to do this, as I’ve run into issues with what should have been simple approaches. An example image is below, and can be downloaded here: https://gofile.io/d/YDAs0G

enter image description here

I’ve tried using both raster’s polygonize tools & stars contouring to no avail. I use a common dataset for all approaches and then outline them below:

s_r <- raster("test_shoreline.tif")

Approach 1:
Use stars st_contour to get the contour of the image. Issue: result is disconnected line segments especially on sharp curves. The missing CRS in sr_lines is not an issue, can reimpose it.

  s_r_st <- st_as_stars(s_r)
  s_r_lines <- st_contour(s_r_st, contour_lines = T, breaks = c(0, 1, 2))
> sr_lines
Simple feature collection with 124 features and 1 field
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 0 ymin: 0.2547393 xmax: 1 ymax: 0.4751185
CRS:            unknown
First 10 features:
   layer                       geometry
1  [0,1] LINESTRING (1 0.4751185, 0....
2  [0,1] LINESTRING (0 0.4751185, 0....
3  [0,1] LINESTRING (0.1196682 0.472...
4  [0,1] LINESTRING (0.1978673 0.472...
5  [0,1] LINESTRING (0.7902844 0.472...
6  [0,1] LINESTRING (0.721564 0.4727...
7  [0,1] LINESTRING (0.7191943 0.468...
8  [0,1] LINESTRING (0.7168246 0.465...
9  [0,1] LINESTRING (0.3613744 0.463...
10 [0,1] LINESTRING (0.7120853 0.463...
plot(sr_lines)

enter image description here

Approach 2:
Use raster rasterToContour to get a linestring. No go as full of multilinestrings.

  s_r_c <- rasterToContour(s_r, nlevels = 2)
  s_r_st <- st_as_sf(s_r_c) 
> s_r_st
Simple feature collection with 1 feature and 1 field
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: 0.001184834 ymin: 0.254737 xmax: 0.9988152 ymax: 0.4751209
CRS:            unknown
    level                       geometry
C_1     1 MULTILINESTRING ((0.0011848...
> st_cast(s_r_st, "LINESTRING")
Simple feature collection with 230 features and 1 field
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 0.001184834 ymin: 0.254737 xmax: 0.9988152 ymax: 0.4751209
CRS:            unknown
First 10 features:
      level                       geometry
C_1       1 LINESTRING (0.001184834 0.4...
C_1.1     1 LINESTRING (0.08412322 0.47...
C_1.2     1 LINESTRING (0.08649052 0.47...
C_1.3     1 LINESTRING (0.1196682 0.472...
C_1.4     1 LINESTRING (0.1220355 0.475...
C_1.5     1 LINESTRING (0.1978673 0.472...
C_1.6     1 LINESTRING (0.2002346 0.475...
C_1.7     1 LINESTRING (0.3234597 0.404...
C_1.8     1 LINESTRING (0.3258294 0.401...
C_1.9     1 LINESTRING (0.325827 0.4063...

enter image description here

Approach 3: Turn the pixels into points and then convert the points into a line. This seems like it should work but the ordering of the coordinates is off?

 s_r_obj <- stars::st_as_stars(s_r) %>%
    st_as_sf(as_points = T, connect8 = T)
  
  s_r_obj2 <- s_r_obj %>%
    dplyr::group_by(layer) %>%
    dplyr::summarise(do_union=T)
  s_r_obj2 <- s_r_obj2[s_r_obj2$layer==1, ]

    s_r_obj3 <- st_cast(s_r_obj2, "LINESTRING")

> s_r_obj3
Simple feature collection with 1 feature and 1 field
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 0.001184834 ymin: 0.2547393 xmax: 0.9988152 ymax: 0.4751185
CRS:            unknown
# A tibble: 1 x 2
  layer                                                                          geometry
  <dbl>                                                                  <LINESTRING [°]>
1     1 (0.001184834 0.4751185, 0.003554502 0.4751185, 0.005924171 0.4751185, 0.00829383~

enter image description here

I was thinking to turn this into a single polygon and then skeletonize as implemented here, but I get multiple polygons from a modified approach 3. Is it simpler to skip the GIS functionality and do a guided-walk along the pixels?

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP