2

I have a number of polygons in 3d from a geojson file, and I would like to make an elevation model. This means that I want a raster, where every pixel is the height of the polygon in this position.

I tried looking at gdal_rasterize, but the description says

As of now, only points and lines are drawn in 3D.

gdal_rasterize

Little geek
  • 592
  • 8
  • 22
  • Can you provide an example of the GeoJSON in the original question? Is each polygon 2D planar but at a different Z value? – Logan Byers Nov 03 '16 at 13:57

1 Answers1

2

I ended up using the scipy.interpolat-function called griddata. This uses a meshgrid to get the coordinates in the grid, and I had to tile it up because of memory restrictions of meshgrid.

import scipy.interpolate as il #for griddata
# meshgrid of coords in this tile
gridX, gridY = np.meshgrid(xi[c*tcols:(c+1)*tcols], yi[r*trows:(r+1)*trows][::-1])

## Creating the DEM in this tile
zi = il.griddata((coordsT[0], coordsT[1]), coordsT[2], (gridX, gridY),method='linear',fill_value = nodata) # fill_value to prevent NaN at polygon outline

The linear interpolation seems to do exactly what I want. See description at https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html

Little geek
  • 592
  • 8
  • 22
  • Would you be willing to share your code starting from the 3d polygon input and resulting in the floating point raster output? Thanks! – Brent Edwards Oct 22 '18 at 20:38