terra dropping / not accepting Z coordinates

[dimfalk]

Some time ago there was an interesting discussion on the r-sig-geo mailing list (Vol. 227, Issue 7+8) about height transformation where @rsbivand pointed out that terra drops resp. does not accept Z coordinates so that e.g. existing POINT Z objects are reduced to POINT objects when making use of terra.

I was waiting for the original questioner to address this here but it seems like no issue was submitted so far. Since my impression was that terra::vect() converts sf objects without any loss of information, I thought I'd mention this issue.

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
library(terra)
#> terra 1.6.17

data <- "Name,Longitude,Latitude,Ellipsoidal_height
gcp1,2.32152788,41.75758799,1235.577
gcp2,2.32129905,41.75757489,1234.38
gcp3,2.3212163,41.75760007,1233.49
gcp4,2.32165884,41.75765011,1235.785
gcp5,2.32163868,41.75764588,1235.741
gcp6,2.32167394,41.75793979,1236.763
gcp7,2.32118573,41.75769993,1232.364
gcp8,2.32115307,41.75770638,1231.897
gcp9,2.32097163,41.75779261,1228.454"

data <- read.table(text = data, header = TRUE, sep = ",")

# sf object representation
data_sf <- st_as_sf(data, 
                    coords = c("Longitude", "Latitude", "Ellipsoidal_height"))

st_geometry(data_sf)
#> Geometry set for 9 features 
#> Geometry type: POINT
#> Dimension:     XYZ
#> Bounding box:  xmin: 2.320972 ymin: 41.75757 xmax: 2.321674 ymax: 41.75794
#> CRS:           NA
#> First 5 geometries:
#> POINT Z (2.321528 41.75759 1235.577)
#> POINT Z (2.321299 41.75757 1234.38)
#> POINT Z (2.321216 41.7576 1233.49)
#> POINT Z (2.321659 41.75765 1235.785)
#> POINT Z (2.321639 41.75765 1235.741)

# SpatVector object representation
data_terra <- vect(data_sf)
geom(data_terra, wkt = TRUE)
#> [1] "POINT (2.321528 41.757588)" "POINT (2.321299 41.757575)"
#> [3] "POINT (2.321216 41.7576)"   "POINT (2.321659 41.75765)" 
#> [5] "POINT (2.321639 41.757646)" "POINT (2.321674 41.75794)" 
#> [7] "POINT (2.321186 41.7577)"   "POINT (2.321153 41.757706)"
#> [9] "POINT (2.320972 41.757793)"

vect(data, geom = c("Longitude", "Latitude", "Ellipsoidal_height"))
#> Error: [vect] the length of 'geom' must be 1 or 2

Would it make sense to extend the existing coordinate representation for Z coordinates according to OGC specifications for WKT? Or at least to throw a warning when Z coordinates are dropped using vect()?

Thank you very much in advance!

[rsbivand]

The R-sig-geo thread is at https://stat.ethz.ch/pipermail/r-sig-geo/2022-July/029040.html

[dimfalk]

Issues seem also to exist with SpatRaster objects with some VERTCRS definition (c.f. this SO question).

Probably there is a way to handle this better than I did, but I feel like this is also connected to this issue about how to properly handle Z coordinates when working with spatial data.

This geodetic perspective focussing on the Z dimension is also partly new to me and I'm not really sure about this, but maybe it's time to think of spatial data in 3D per default? I feel like focus in the past (completely subjective) was almost entirely on 2D and this might be changing recently...

[rhijmans]

Thanks, I do not think these two issues ("SpatVector ignoring z-values" and "crs transformation error") are in any way related

[dimfalk]

Thanks for your reply, Robert!

I don't think these two issues to be caused by the same underlying problem either. I was just wondering if the solution for both might not be a similar one - to extend the spatial model in z dimension - hence the mention.

Would the "crs transformation error" part be worth a separate issue from your point of view?

[rsbivand]

@dimfalk 4326 appears to say that it is "Horizontal component of 3D system.":

$ projinfo EPSG:4326
PROJ.4 string:
+proj=longlat +datum=WGS84 +no_defs +type=crs

WKT2:2019 string:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

But maybe EPSG:4979:

$ projinfo EPSG:4979
PROJ.4 string:
+proj=longlat +datum=WGS84 +no_defs +type=crs

WKT2:2019 string:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,3],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["ellipsoidal height (h)",up,
            ORDER[3],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Geodesy. Navigation and positioning using GPS satellite system."],
        AREA["World: Afghanistan, Albania, Algeria, American Samoa, Andorra, Angola, Anguilla, Antarctica, Antigua and Barbuda, Argentina, Armenia, Aruba, Australia, Austria, Azerbaijan, Bahamas, Bahrain, Bangladesh, Barbados, Belgium, Belgium, Belize, Benin, Bermuda, Bhutan, Bolivia, Bonaire, Saint Eustasius and Saba, Bosnia and Herzegovina, Botswana, Bouvet Island, Brazil, British Indian Ocean Territory, British Virgin Islands, Brunei Darussalam, Bulgaria, Burkina Faso, Burundi, Cambodia, Cameroon, Canada, Cape Verde, Cayman Islands, Central African Republic, Chad, Chile, China, Christmas Island, Cocos (Keeling) Islands, Comoros, Congo, Cook Islands, Costa Rica, Côte d'Ivoire (Ivory Coast), Croatia, Cuba, Curacao, Cyprus, Czechia, Denmark, Djibouti, Dominica, Dominican Republic, East Timor, Ecuador, Egypt, El Salvador, Equatorial Guinea, Eritrea, Estonia, Eswatini (Swaziland), Ethiopia, Falkland Islands (Malvinas), Faroe Islands, Fiji, Finland, France, French Guiana, French Polynesia, French Southern Territories, Gabon, Gambia, Georgia, Germany, Ghana, Gibraltar, Greece, Greenland, Grenada, Guadeloupe, Guam, Guatemala, Guinea, Guinea-Bissau, Guyana, Haiti, Heard Island and McDonald Islands, Holy See (Vatican City State), Honduras, China - Hong Kong, Hungary, Iceland, India, Indonesia, Islamic Republic of Iran, Iraq, Ireland, Israel, Italy, Jamaica, Japan, Jordan, Kazakhstan, Kenya, Kiribati, Democratic People's Republic of Korea (North Korea), Republic of Korea (South Korea), Kosovo, Kuwait, Kyrgyzstan, Lao People's Democratic Republic (Laos), Latvia, Lebanon, Lesotho, Liberia, Libyan Arab Jamahiriya, Liechtenstein, Lithuania, Luxembourg, China - Macao, Madagascar, Malawi, Malaysia, Maldives, Mali, Malta, Marshall Islands, Martinique, Mauritania, Mauritius, Mayotte, Mexico, Federated States of Micronesia, Monaco, Mongolia, Montenegro, Montserrat, Morocco, Mozambique, Myanmar (Burma), Namibia, Nauru, Nepal, Netherlands, New Caledonia, New Zealand, Nicaragua, Niger, Nigeria, Niue, Norfolk Island, North Macedonia, Northern Mariana Islands, Norway, Oman, Pakistan, Palau, Panama, Papua New Guinea (PNG), Paraguay, Peru, Philippines, Pitcairn, Poland, Portugal, Puerto Rico, Qatar, Reunion, Romania, Russian Federation, Rwanda, St Barthelemy, St Kitts and Nevis, St Helena, Ascension and Tristan da Cunha, St Lucia, St Martin, St Pierre and Miquelon, Saint Vincent and the Grenadines, Samoa, San Marino, Sao Tome and Principe, Saudi Arabia, Senegal, Serbia, Seychelles, Sierra Leone, Singapore, Slovakia (Slovak Republic), Slovenia, St Maarten, Solomon Islands, Somalia, South Africa, South Georgia and the South Sandwich Islands, South Sudan, Spain, Sri Lanka, Sudan, Suriname, Svalbard and Jan Mayen, Sweden, Switzerland, Syrian Arab Republic, Taiwan, Tajikistan, United Republic of Tanzania, Thailand, The Democratic Republic of the Congo (Zaire), Togo, Tokelau, Tonga, Trinidad and Tobago, Tunisia, Türkiye (Turkey), Turkmenistan, Turks and Caicos Islands, Tuvalu, Uganda, Ukraine, United Arab Emirates (UAE), United Kingdom (UK), United States (USA), United States Minor Outlying Islands, Uruguay, Uzbekistan, Vanuatu, Venezuela, Vietnam, US Virgin Islands, Wallis and Futuna, Western Sahara, Yemen, Zambia, Zimbabwe."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4979]]

which has a vertical axis? The ellipsoidal height may not be appropriate for the case, though.

[rsbivand]

I get (settng PROJ_NETWORK=ON before starting R):

> dsm
class       : SpatRaster 
dimensions  : 1500, 1501, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 25349, 26850, 170258, 171758  (xmin, xmax, ymin, ymax)
coord. ref. : NAD83(2011) / WISCRS Brown (m) (EPSG:7532) 
source(s)   : memory
name        :       Z 
min value   : 185.836 
max value   : 333.709 
> library(sf)
Linking to GEOS 3.11.1, GDAL 3.6.0, PROJ 9.1.1; sf_use_s2() is TRUE
> sf_proj_network()
[1] TRUE
> sf_proj_pipelines(crs(dsm), "EPSG:4979")
Candidate coordinate operations found:  5 
Strict containment:     FALSE 
Axis order auth compl:  FALSE 
Source:  COMPOUNDCRS["NAD83(2011) / WISCRS Brown (m) + NAVD88 height - Geoid18 (m)",
    PROJCRS["NAD83(2011) / WISCRS Brown (m)",
        BASEGEOGCRS["NAD83(2011)",
            DATUM["NAD83 (National Spatial Reference System 2011)",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",6318]],
        CONVERSION["unnamed",
            METHOD["Transverse Mercator",
                ID["EPSG",9807]],
            PARAMETER["Latitude of natural origin",43,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",-88,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["Scale factor at natural origin",1.00002,
                SCALEUNIT["unity",1],
                ID["EPSG",8805]],
            PARAMETER["False easting",31600,
                LENGTHUNIT["meter",1],
                ID["EPSG",8806]],
            PARAMETER["False northing",4600,
                LENGTHUNIT["meter",1],
                ID["EPSG",8807]]],
        CS[Cartesian,2],
            AXIS["x",east,
                ORDER[1],
                LENGTHUNIT["meter",1]],
            AXIS["y",north,
                ORDER[2],
                LENGTHUNIT["meter",1]],
        ID["EPSG",7532]],
    VERTCRS["NAVD88 height - Geoid18 (m)",
        VDATUM["North American Vertical Datum 1988"],
        CS[vertical,1],
            AXIS["up",up,
                LENGTHUNIT["meter",1]],
        GEOIDMODEL["GEOID18"],
        ID["EPSG",5703]]] 
Target:  EPSG:4979 
Best instantiable operation has accuracy: 2.015 m
Description: Inverse of unnamed + Inverse of NAD83(2011) to NAVD88 height
             (3) + NAD83(2011) to WGS 84 (1) + axis order
             change (2D)
Definition:  +proj=pipeline +step +inv +proj=tmerc +lat_0=43 +lon_0=-88
             +k=1.00002 +x_0=31600 +y_0=4600 +ellps=GRS80 +step
             +proj=vgridshift +grids=us_noaa_g2018u0.tif
             +multiplier=1 +step +proj=unitconvert +xy_in=rad
             +xy_out=deg
> dsm_test <- terra::project(dsm, "EPSG:4979", method="bilinear")
> dsm_test
class       : SpatRaster 
dimensions  : 1235, 1727, 1  (nrow, ncol, nlyr)
resolution  : 1.093727e-05, 1.093727e-05  (x, y)
extent      : -88.0786, -88.05972, 44.49092, 44.50443  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4979) 
source(s)   : memory
name        :        Z 
min value   : 185.8360 
max value   : 295.1357 

but I'm afraid that Z is not handled at all.

[rhijmans]

@dimfalk, yes this would have been good for a new issue.

[dimfalk]

@rsbivand: Wow, I feel like there is a lot more to learn about this. Did not knew about EPSG:4979 until now, also thanks for hinting at sf_proj_pipelines() - looks promising! But I assume this discussion is better to be continued in the new issue since these issues are not related to keep thing on-topic here?

@rhijmans: Will do!

[twest820]

In #125 back in January 2021, @rhijmans wrote

Thanks, this example now works; although terra indeed ignores Z/M values but that should change soon.

which in February 2021 changed to

It would, of course, be better to just read the Z/M values. And I should eventually do that. But I am not too inspired as I never use them, so if anyone out there would really like to have them, let me know.

As of today, with terra 1.7.29 if I do reducedPointCloud = vect("marked 3D points in xyz format.gpkg") I get

[vect] Z coordinates ignored 

and geom(derivedPointCloud ) predictably has only x and y coordinates. This means I've literally gigabytes of data that can't be usefully loaded using vect() because it's written out in standard xyz format rather than pushing the points' z coordinates into the attribute table.

I guess this puts me in the category of people who'd really like to have z support, hence the bump on this open issue. Pretty much all of the point processing I've done in the past several months is 3D and a lot of it's actually marked point processes. So I've a number of other use cases where wkbPointMZ or wkbPointM formats make sense alongside wkbGeometryType.wkbPoint25D—only reason I'm not using measure formats is it's helpful to make the files self documenting by having an attribute table column containing the name of what measurement the points are marked with but z is standard enough that's not a concern.

[andyatsenversa]

This isn't a solution, but a workaround that has worked for my current scenario (lots of dwgs with 3d lines i need to extract). TLDR: read the geometry as WKT using OGR-SQL, then use that to rebuild the original xyz geometry.

The workaround is possible due to terra::vect accepting OGR-SQL dialect options for handling geometry, including the ability to select geometry as Well-Known-Text (OGR_GEOM_WKT): https://gdal.org/user/ogr_sql_dialect.html

library(dplyr)
library(sf)
library(terra)

# define the project coordinate system
projcrs <- 4326

# define a path to the source dwg file
file_src <- "path_to_survey.dwg"

# extract the list of layers in the source file
map_layers <- terra::vector_layers(file_src)

# iterate through each layer, build a query that selects all fields (*) and also 'OGR_GEOM_WKT' 
extract_dwg_3d <- lapply(map_layers, function(y) {
        vect(file_src, query = paste0("SELECT OGR_GEOM_WKT as geometry, * from ", '"', y, '"')) |> # read one layer
            as_tibble() |> # drop the terra geometry
            mutate(layer = y) |> # add the layer name
            st_as_sf(wkt = "geometry", crs = projcrs) # use sf to rebuild the geometry based on the WKT with XYZ
    }) |> bind_rows()