3

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 below 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? The image below effectively illustrates what I am trying to achieve. 

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)

enter image description here

Improve this question
5
  • Nice question and I wish everyone new to GIS had problems like this. Could you just explain what you mean by "relative scale is lost"?
    – John Powell
    Commented Nov 8, 2016 at 14:50
  • Hi John, here's what I suspect the problem is: the raster (as seen on left) upper left hand's coordinates are (x1, y1) whereas the building is at (x2, y2). When I get the geom from PostGIS as a raster, I try to load it into a gdal driver where I have set the geo_transform and projection of the main raster.
    – sailor
    Commented Nov 8, 2016 at 18:29
  • I would hope that the building raster would then be projected into the space of the main raster, so we would get a raster starting at (x1, y1) that would be empty except for building starting at (x2, y2). However, what really happens (visually, as seen in the display) is that I get a building that takes up the whole panel above. Nor do I know if it's actually at (x1, y1).
    – sailor
    Commented Nov 8, 2016 at 18:29
  • If I understand well, you want one raster for each building. Each one of these rasters have to be the size (width and height) of the building raster (on the left), contain only one building each and the buildings have to be alingneg with the corresponding ones in the building raster?
    – Pierre Racine
    Commented Nov 8, 2016 at 22:48
  • Hi Pierre, that's close. Everything you have said is accurate, except I would add that several buildings could be a single raster (as seen on right). Ultimately, the goal is to get a numpy array with 1's where there all buildings and 0's where there are not.
    – sailor
    Commented Nov 9, 2016 at 1:50

1 Answer 1

Reset to default
1

ST_AsRaster() produce one raster per polygon. If you want all the buidings in the same raster at the end, you have to ST_Union(rast, 'max') all of them together. Try this:

SELECT ST_AsGDALRaster(ST_Union(ST_AsRaster(geom, 10, 10, '16BUI'), 'max'), GTiff) 
 FROM buildings 
WHERE ST_ContainsProperly(
          ST_Polygon(
              (SELECT rast 
                 FROM eightbands WHERE rid = %s), 1), 
          geom);

Please refer to this guide for proper and efficient rasterization of a geometry table.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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