I am new to GIS and am trying to build an classifier for buildings using labeled satellite imagery.

I have a PostGIS database where I am storing rasters and POLYGON geom's which represent buildings. My goal is to create a numpy array for each raster where buildings are represented as 1's and everything else as 0's.

I am able to make this work by running pgsql2shp on each geom; opening a gdal driver and setting the geo_transform and projection on the driver to the same as the raster's; and then rasterizing the layer (burning the values). However, this involves an intermediate step of outputing a shapefile. I would prefer to extract each geom directly from PostGIS using a query and project onto the raster, but it seems that the relative scale is lost when I use ST_AsRaster. See [A] for the query I was using (I tried a few variants of ST_AsRaster, this is just to illustrate my general approach).

How can get out a properly scaled raster (or some other format I can convert to a numpy array) from my database? [B] effectively illustrates what I trying to achieve.

[A] "SELECT ST_AsGDALRaster(ST_AsRaster(geom, 10, 10, '16BUI'), \

•    'GTiff') FROM buildings WHERE ST_ContainsProperly \
  

•    (ST_Polygon((SELECT rast FROM eightbands WHERE rid = %s), 1), \
  

•    geom)"
  

[B]