Add Z values to preexisting geometry

[seanreilly66]

Hi,

Thank you for your work on this package. I have a relatively simple question with which I would appreciate your help.

I have an xy point sf object to which I am hoping to add height values from an external vector. I feel like this should be an easy operation but I have not been able to figure out the best way to do so. Any advice would be greatly appreciated.

Here is a demonstration:

library(sf)

sf1 <- data.frame(x = runif(4), y = runif(4)) |>
  st_as_sf(coords = c('x', 'y'))

external_z <- runif(4)

sf1$z <- external_z

# Update geometry in sf1 to xyz with the Z values from "external_z"

[twest820]

Hit this question too. The answer I was able to find is more complicated, and thus substantially slower, than copying elements of vector into an array of structs should be. But at least it works.

library(dplyr)

sf1 = st_zm(sf1, drop = FALSE, what = "Z")
sf1xyz = st_coordinates(sf1)
sf1xyz[, "Z"] = external_z
sf1xyz = as_tibble(sf1xyz) %>% rowwise() %>% mutate(geom = list(st_point(c(X, Y, Z))))
sf1xyz = st_sfc(sf1xyz$geom)
sf1 = st_set_geometry(sf1, sf1xyz)

Looks to me like there's three basic feature gaps introducing complexity and performance issues here.

1. st_zm(drop = FALSE) doesn't accept a z (or m) vector to add.
2. st_set_geometry() doesn't accept matrices even though they're what st_coordinates() returns. Or, equivalently, there is no st_set_coordinates().
3. st_as_sfc() doesn't accept matrices either, even though its purpose is to do what's needed here.

A corollary issue is I didn't have luck finding a lighter weight way than dplyr's rowwise() to yield a list out of calling a function on each row of a matrix. R's apply family of functions seems to lack a variant for this, instead doing weird stuff like silently flattening st_points back to their coordinates. Maybe there's something less cumbersome in data.table.

The way I encountered this in real world coding was

merge = bind_rows(xyLayerWithCompoundCrs,
                  st_transform(xyzLayer1withProjectedCrs, crs(xyLayerWithCompoundCrs)),
                  st_transform(xyzLayer2withProjectedCrs, crs(xyLayerWithCompoundCrs)))
merge
Simple feature collection with 77049 features and 10 fields
Geometry type: POINT
Dimension:     XY, XYZ

st_coordinates(merge)
Error in dimnames(x) <- dn : 
  length of 'dimnames' [2] not equal to array extent
In addition: Warning message:
In matrix(unlist(x, use.names = FALSE), nrow = length(x), byrow = TRUE,  :
  data length [211178] is not a sub-multiple or multiple of the number of rows [77049]

Which illustrates a couple related issues,

4. st_transform() doesn't transform xy to xyz points when transforming from a 2D to a 3D CRS. Seems like there should be at least a warning about this.
5. bind_rows() doesn't make its output consistent. There should probably be a warning on this like there is when binding CRSes that don't match.

Also, it's kinda odd it seems mandatory to store points' coordinates in an array of structs (AoS) instead of a structure of arrays (SoA). Usually SoA is preferred in new code for fewer memory allocations and SIMD friendliness. For example, AoS requires allocating hundreds of thousands of temporary point objects here just to copy three scalars into and out of each of them once.

[edzer]

The way to do this is to have the Z coordinate in the data.frame and then convert that into an sf object, e.g.

library(sf)
# Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
data.frame(x = 1:3, y = 1:3, z = 1:3, a = rnorm(3)) |>
  st_as_sf(coords = c("x", "y", "z"))
# Simple feature collection with 3 features and 1 field
# Geometry type: POINT
# Dimension:     XYZ
# Bounding box:  xmin: 1 ymin: 1 xmax: 3 ymax: 3
# CRS:           NA
#            a        geometry
# 1 -0.7473629 POINT Z (1 1 1)
# 2 -0.1859593 POINT Z (2 2 2)
# 3 -0.1143853 POINT Z (3 3 3)

@twest820 thanks for your comments! On your point 4: does it work if the coordinates are POINT Z geometries with zero Z values? on your point 5: bind_rows() is not a generic, so sf cannot do anything to make this work. On your last point: this is well known, and a consequence of having sf POINT geometries being implemented similar to other geometries, which is quite convenient. See also #2376 and the pointx branch (which is unlikely to be merged any time soon).

[twest820]

st_as_sf() unfortunately isn't viable for adding z values to a preexisting geometry because st_set_geometry() doesn't support sf.

sf1 = st_zm(sf1, drop = FALSE, what = "Z")
sf1 = st_set_geometry(sf1, st_as_sf(as_tibble(st_coordinates(sf1)) %>% mutate(Z = external_z), coords = c("X", "Y", "Z")))
Error in `st_geometry<-.sf`(`*tmp*`, value = value) : 
  is.character(value) || inherits(value, "sfc") is not TRUE

Per that error, workable data flow with sf 1.0-19 seems to be restricted to st_sfc().

Yes, bind_rows() does fine if its inputs are made minimally consistent.

merge = bind_rows(st_zm(xyLayerWithCompoundCrs, drop = FALSE, what = "Z"),
                  st_transform(xyzLayer1withProjectedCrs, crs(xyLayerWithCompoundCrs)),
                  st_transform(xyzLayer2withProjectedCrs, crs(xyLayerWithCompoundCrs)))

[edzer]

You need to take the st_geometry() from the output of st_as_sf().

[twest820]

Ah, thanks. If it's of interest here's the basic flow for my real world use.

merge1 = bind_rows(...)
merge2 = bind_rows(st_transform(st_zm(...), crs = crs(merge1), ...)
merge2xyz = st_geometry(merge2)
...
output = st_set_geometry(merge1[merge1indices, ], st_geometry(st_as_sf(merge2xyz[merge2indices, c("X", "Y")], coords = c("X", "Y"), crs = crs(merge1))))

Other readers may find it helpful to note st_set_geometry() is really st_set_geometry_and_crs() since it slices the CRS off the merge1 slice if a CRS isn't passed to st_as_sf().