[question] M and XYZM are suported for `sfg` but not `sfc`

[JosiahParry]

Hi! I am in the process of writing conversions to/from EsriJSON geometry objects for sfg, sfc, and sf objects. I've completed the sfg conversions. However, I am now noticing that M and ZM are not supported in sfc objects. Is there any documentation on why this is or how to get around handling it?

(mpnt <- sf::st_multipoint(matrix(rnorm(12), 3)))
#> MULTIPOINT ZM ((1.361224 -1.337553 -0.4722005 1.366485), (1.724852 1.249042 0.525239 0.7809072), (1.200532 0.1311064 0.2935638 -0.487114))

sf::st_sfc(mpnt) |>
  sf::st_cast("POINT")
#> Error in eval(expr, envir, enclos): GEOS does not support XYM or XYZM geometries; use st_zm() to drop M

^{Created on 2024-04-29 with reprex v2.1.0}

[edzer]

They are supported, up to casting them:

library(sf)
# Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
(mpnt <- sf::st_multipoint(matrix(rnorm(12), 3)))
# MULTIPOINT ZM ((1.261043 -0.08883676 0.5119583 2.707271), (0.3177616 -0.2379901 0.02654004 0.5705929), (0.5544696 0.5938268 -1.428157 -0.1402494))
sf::st_sfc(mpnt)
# Geometry set for 1 feature 
# Geometry type: MULTIPOINT
# Dimension:     XYZM
# Bounding box:  xmin: 0.3177616 ymin: -0.2379901 xmax: 1.261043 ymax: 0.5938268
# z_range:       zmin: -1.428157 zmax: 0.5119583
# m_range:       mmin: -0.1402494 mmax: 2.707271
# CRS:           NA
# MULTIPOINT ZM ((1.261043 -0.08883676 0.5119583 ...

the annoying thing is that st_cast() calls st_is_empty(), which calls the GEOS function:

4: CPL_geos_is_empty(st_geometry(x))
3: st_is_empty(x)
2: st_cast.sfc(x, "POINT")
1: st_cast(x, "POINT")

which doesn't like ZM: Z is not the problem:

library(sf)
# Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
(mpnt <- sf::st_multipoint(matrix(rnorm(12), 4)))
# MULTIPOINT Z ((0.5337231 -0.3189474 0.878158), (0.2367451 2.446873 -0.7087885), (0.2484128 1.026354 1.515804), (-0.7179272 1.957629 -0.3093828))
sf::st_sfc(mpnt) |> st_cast("POINT")
# Geometry set for 4 features 
# Geometry type: POINT
# Dimension:     XYZ
# Bounding box:  xmin: -0.7179272 ymin: -0.3189474 xmax: 0.5337231 ymax: 2.446873
# CRS:           NA
# POINT Z (0.5337231 -0.3189474 0.878158)
# POINT Z (0.2367451 2.446873 -0.7087885)
# POINT Z (0.2484128 1.026354 1.515804)
# POINT Z (-0.7179272 1.957629 -0.3093828)

but it dislikes M. One could program around the st_is_empty(), but st_cast() already is a big mess as it stands now.

[JosiahParry]

Oh, interesting. The following works as a work around:

matrix(rnorm(12), 3) |> 
  apply(1, sf::st_point, simplify = FALSE) |> 
  sf::st_sfc()
#> Geometry set for 3 features 
#> Geometry type: POINT
#> Dimension:     XYZM
#> Bounding box:  xmin: -0.933232 ymin: 0.498481 xmax: 2.199169 ymax: 1.742439
#> z_range:       zmin: -1.95358 zmax: 0.9565412
#> m_range:       mmin: -1.878737 mmax: -0.02329908
#> CRS:           NA
#> POINT ZM (1.138996 1.742439 0.9565412 -0.02329908)
#> POINT ZM (2.199169 0.8638457 -1.95358 -1.878737)
#> POINT ZM (-0.933232 0.498481 -0.08867432 -1.127...

^{Created on 2024-04-29 with reprex v2.1.0}

[edzer]

Yes, but this pipeline doesn't work:

st_multipoint() |> 
  apply(1, sf::st_point, simplify = FALSE) |> 
  sf::st_sfc()

which is the point in case you stumbled on.