I have Open Street Map [OSM] data imported so that I can access the geometry of buildings. I have also loaded SRTM digital height data from NASA from a GeoTIFF raster. I have projected the raster into EPSG:3857 (the same as OSM) from the original EPSG:4326.

For the SRTM tiles that I am interested in, I have warped them into 1 meter resolution with cubic spline interpolation via ST_Resample and produced smaller tiles from that with ST_Tile(as a single INSERT .. SELECT which is very handy).

Now I run into a problem:

The ST_Intersection() of the resampled tile and the building polygon simply produces a square bounding box, which is not what I want. I want to get at the exact x,y points on the raster for the intersection because I'm using this data for generating 3d voxel scenes (1 pixel = one voxel) and I want to produce a 2d map of the rendered area to pixel-point precision as well.

What I am currently doing:

I use ST_Clip() to clip the raster by the geometry of the buildings, which puts NODATA data into the returned clipped raster everywhere except where the building is actually situated. This is like using an template to "burn" the image of the building into the raster, leaving the points inside the building completely intact while everything else is burned away.

I dump the resulting clipped raster with ST_PixelAsPoints() or ST_PixelAsCentroids() which gives me x,y positions on the raster where the building is located.

Is this the right approach?

I have heard of something called a 'pointcloud' but PostGIS doesn't seem to have this available. Converting the raster into polygons doesn't work, because the x,y data from the raster is lost.

  • 1
    Welcome to gis.SE. You are using SRTM, which is around 90m horizontal per sample. I don't think you're going to get "pixel precision" at any meaningful building scale out of that data. – BradHards Jan 25 '14 at 22:52
  • I interpolated it to a higher resolution. This is like zooming into a picture. The results are not 100% accurate but they look pretty good. The same concept applies to any size region. I could be clipping a county to the (unscaled) raster, so I don't think your suggestion applies. What I'm doing is works, and produces reasonable results, but is rather computationally intensive and I was wondering if there was a better way. – Justin Swanhart Jan 26 '14 at 5:38
  • Note, I'm also determining on which pixels roads and other items overlap so that I can produce the Z order properly. I don't think SQLFiddle has 9.3 for an example. Here is a a clipped raster so you can see what I mean by "point precision" of buildings near the rim of the grand canyon: picpaste.com/sample_elevation_map.png – Justin Swanhart Jan 26 '14 at 6:48

Your Answer

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

Browse other questions tagged or ask your own question.