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.