1

I am trying to convert some DEM in .tif image to polygon geometries representing 3 dimensional triangles that is Well-known text (WKT) - 3D (x, y, z)

I was able to use SAGA to convert my GeoTIFF image to TIN, but what really wanted is a WKT format as stated above.

How do I do this is in QGIS, SAGA, GDAL, PostGIS or any other open source solution?

|improve this question
  • If you've already done the first step (GeoTIFF -> TIN), then all you really need is TIN -> WKT, which at least some, probably all, of those packages can accomplish. Please choose one package and attempt the conversion, then if you have a problem, edit the question with details of that problem. – Vince Nov 9 '18 at 9:56
1

In R you can do this with the dev-version of quadmesh, and the sf package. 

library(raster)
## create a raster 
r <- raster(volcano)

## convert to indexed quad form, the rgl mesh3d format
qm <- quadmesh::quadmesh(r)

## split quad index into triangle index
tri <- quadmesh::triangulate_quads(qm$ib)

## convert to simple features
ps <- lapply(split(rbind(tri, tri[1, ]), rep(seq_len(ncol(tri)), each = 4)), 
       function(idx) sf::st_polygon(list(t(qm$vb[1:3, idx]))))
p_xyz <- sf::st_sf(geometry = sf::st_sfc(ps, crs = projection(r)))

p_xyz

Simple feature collection with 10614 features and 0 fields
geometry type:  POLYGON
dimension:      XYZ
bbox:           xmin: 1.490116e-08 ymin: 1.490116e-08 xmax: 1 ymax: 1
epsg (SRID):    NA
proj4string:    NA
First 10 features:
                         geometry
1  POLYGON Z ((1.490116e-08 1 ...
2  POLYGON Z ((0.01639344 1 99...
3  POLYGON Z ((0.01639344 1 99...
4  POLYGON Z ((0.03278689 1 10...
5  POLYGON Z ((0.03278689 1 10...
6  POLYGON Z ((0.04918033 1 10...
7  POLYGON Z ((0.04918033 1 10...
8  POLYGON Z ((0.06557377 1 10...
9  POLYGON Z ((0.06557377 1 10...
10 POLYGON Z ((0.08196721 1 10...

Then convert to text form with 

sf::st_as_text(sf::st_geometry(p_xyz))

For triangulate_quads the development version is required, e.g. 

devtools::install_github("hypertidy/quadmesh")

https://github.com/hypertidy/quadmesh

|improve this answer
  • I don't much of R. But I will try it out. – Umar Yusuf Nov 9 '18 at 20:42
  • > tri <- quadmesh::triangulate_quad(qm$ib) I couldn't pass the line above for this.... "Error: 'triangulate_quad' is not an exported object from 'namespace:quadmesh'" I not sure what package is causing that error? My sessionInfo() is.... – Umar Yusuf Nov 9 '18 at 22:05
  • R version 3.5.1 (2018-07-02) Platform: x86_64-w64-mingw32/x64 (64-bit) other attached packages: [1] raster_2.8-4 sp_1.3-1 loaded via a namespace (and not attached): [1] Rcpp_0.12.19 lattice_0.20-35 codetools_0.2-15 [4] crayon_1.3.4 dplyr_0.7.7 assertthat_0.2.0 [7] grid_3.5.1 R6_2.3.0 magrittr_1.5 [10] pillar_1.3.0 rlang_0.3.0.1 bindrcpp_0.2.2 [13] tools_3.5.1 glue_1.3.0 purrr_0.2.5 [16] compiler_3.5.1 pkgconfig_2.0.2 quadmesh_0.1.0 [19] tidyselect_0.2.5 bindr_0.1.1 tibble_1.4.2 > – Umar Yusuf Nov 9 '18 at 22:05
  • quadmesh 0.1.0 on CRAN does not have the triangulate_quads function, the dev version is needed, installed from github as per the note above. (I know it's not great having to install dev code, but it's actually so simple it seemed worthwhile - you could just copy the source and use it direct - literally a single line github.com/hypertidy/quadmesh/blob/…). – mdsumner Nov 9 '18 at 23:25
  • Ok, I have succeeded in installing the development version.... I will try it on my local data to see the result. – Umar Yusuf Nov 10 '18 at 15:53
0

If you do not mind mixing SAGA and QGIS operations:

SAGA steps (1, 2, 3)

  1. Grid to TIN or Grid to TIN (Surface Specific Points): convert DEM to TIN (which you have done so already).
  2. TIN to Shapes in (TIN | Tools): exports TIN to Shapefiles. Please select <Points option, which produces 2D Point Shapefile with an attribute Value (height data).
  3. Convert Vertex Type (2D/3D) (in Shapes | Tools): converts 2D Shapefile to 3D Shapefile. Please assign Value field above to the Z option.

QGIS step

  1. Read and Save the Shapefile as a csv, with export geometry as WKT. This will create a 2D WKT field and Value field (separately). Then edit/concatenate these fields by a text editor.
|improve this answer
  • that you for the steps. However, this lead to WKT points like these: POINT ZM (-116.9125692686 32.0029312393 17.7170677185059 0). What I wanted should be the POLYGONS (that is the TIN maybe with three sets of wkt) like: POLYGON Z ((...., ....., ....)? – Umar Yusuf Nov 9 '18 at 20:38

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.