R: convert TIN to multipolygon with package sf

Question

I have a dataset with three-dimensional dataset that has multipolygons and TINs with z-values. I would like to convert these into (multi)polygons.

This is my example data:

element <- st_as_sf( st_as_sfc("TIN(((2661480.805 1241060.614 574.379,2661482.629 1241050.38 574.379,2661482.611 1241050.377 574.379,2661480.805 1241060.614 574.379)),((2661480.785 1241060.611 574.379,2661480.805 1241060.614 574.379,2661482.611 1241050.377 574.379,2661480.785 1241060.611 574.379)),((2661473.445 1241059.303 574.379,2661480.785 1241060.611 574.379,2661475.269 1241049.07 574.379,2661473.445 1241059.303 574.379)),((2661475.269 1241049.07 574.379,2661480.785 1241060.611 574.379,2661482.611 1241050.377 574.379,2661475.269 1241049.07 574.379)),((2661473.414 1241059.298 574.379,2661473.445 1241059.303 574.379,2661475.238 1241049.064 574.379,2661473.414 1241059.298 574.379)),((2661475.238 1241049.064 574.379,2661473.445 1241059.303 574.379,2661475.269 1241049.07 574.379,2661475.238 1241049.064 574.379)))"))

Is there any way to convert this to polygons using R?

Answer 1

Your element is a data frame with the TIN data in element$x:

> element$x
Geometry set for 1 feature 
geometry type:  TIN
dimension:      XY
bbox:           xmin: 2661473 ymin: 1241049 xmax: 2661483 ymax: 1241061
epsg (SRID):    NA
proj4string:    NA
TIN (((2661481 1241061 574.379, 2661483 1241050...

Then element$x is a list with one element in it:

> element$x[[1]]
TIN (((2661481 1241061 574.379, 2661483 1241050 574.379, 2661483 1241050 574.379, 2661481 1241061 574.379)), ((2661481 1241061 574.379, 2661481 1241061 574.379, 2661483 1241050 574.379, 2661481 1241061 574.379)), ((2661473 1241059 574.379, 2661481 1241061 574.379, 2661475 1241049 574.379, 2661473 1241059 574.379)), ((2661475 1241049 574.379, 2661481 1241061 574.379, 2661483 1241050 574.379, 2661475 1241049 574.379)), ((2661473 1241059 574.379, 2661473 1241059 574.379, 2661475 1241049 574.379, 2661473 1241059 574.379)), ((2661475 1241049 574.379, 2661473 1241059 574.379, 2661475 1241049 574.379, 2661475 1241049 574.379)))

which seems to be the triangles of the TIN. You can get each triangle out as list elements:

> element$x[[1]][[1]]
TRIANGLE Z ((2661481 1241061 574.379, 2661483 1241050 574.379, 2661483 1241050 574.379, 2661481 1241061 574.379))
> element$x[[1]][[2]]
TRIANGLE Z ((2661481 1241061 574.379, 2661481 1241061 574.379, 2661483 1241050 574.379, 2661481 1241061 574.379))

And those list elements seem to be 3 column matrices with an extra class. Strip the class off to see the bare matrix:

> unclass(element$x[[1]][[2]])
[[1]]
        [,1]    [,2]    [,3]
[1,] 2661481 1241061 574.379
[2,] 2661481 1241061 574.379
[3,] 2661483 1241050 574.379
[4,] 2661481 1241061 574.379

Those feed nicely into st_polygon:

> st_polygon(unclass(element$x[[1]][[2]]))
POLYGON Z ((2661481 1241061 574.379, 2661481 1241061 574.379, 2661483 1241050 574.379, 2661481 1241061 574.379))

(although you don't need the unclass here or below really)

Wrapping that up in some apply nonsense and doing it all into an SFC vector:

> polys = do.call(st_sfc, lapply(element$x[[1]], function(tin){st_polygon(unclass(tin))}))
> polys
Geometry set for 6 features 
geometry type:  POLYGON
dimension:      XYZ
bbox:           xmin: 2661473 ymin: 1241049 xmax: 2661483 ymax: 1241061
z_range:        zmin: 574.379 zmax: 574.379
epsg (SRID):    NA
proj4string:    NA
First 5 geometries:
POLYGON Z ((2661481 1241061 574.379, 2661483 12...
POLYGON Z ((2661481 1241061 574.379, 2661481 12...
POLYGON Z ((2661473 1241059 574.379, 2661481 12...
POLYGON Z ((2661475 1241049 574.379, 2661481 12...
POLYGON Z ((2661473 1241059 574.379, 2661473 12...

When plotted that only looks like two triangles:

So I guess they overlap? All the Z values are the same too. In a TIN the Z values can be different at each point so you need to do some work to assign a Z to a polygon. Suggest you do st_coordinates(polys) and then work on Z aggregated by L1.

I'm not sure how much this is documented so wrap it all in a function and keep it in one place so if the underlying representation changes you can fix your code easily.