2

I'm having trouble creating a raster as a feature map for my algorithm.

My algorithm needs a numpy array of fixed dimension (e.g. 200x200) indicating if there is land (1) or no land (0) at that position.

I have a PostgreSQL database of Heteogeneous Integrated Dataset https://core.ac.uk/display/158460560 which contains a table named 'geographic_features.europe_coastline_polygone'.

I select geometry for my specified location xmin, ymin, xmax, ymax using a intersection query and this seems to be fine.

But now I want to rasterize the geometry. Since I've found that creating an empty raster and using that as template I'm doing that.

ST_MakeEmptyRaster({width}, {height}, {xmin}, {ymax}, {x_width}, {y_width}, 0, 0, 4326)

All of these parameters are calculated in python. x_width and y_width are calculated as x_width = (xmax - xmin) / width.

Then I use ST_asRaster(geom, template_raster)

The loading is done via WKB using the second answer from here Downloading raster data into python from postgis using psycopg2

My problem is that I'm not getting a array of size (width, height) = (20, 20) as expected, but an array of size (12, 12).

I'm assuming that PostGIS only rasters the bounding box of the geometry instead of the specified width and height. Am I wrong here? Searching the problem only yielded questions about fixing the cell width, which seems to be working fine already for me.

0

First you should always have a visual look at the results of your queries. I use OpenJump which has a very simple interface with PostGIS. Unfortunately it does not support postgis rasters directly but you can always vectorise raster tiles this way:

SELECT rast::geometry, rid
FROM yourRasterTable;

and then view groups of pixel values for one tile this way:

SELECT (gv).geom, (gv).val
FROM (SELECT ST_DumpAsPolygons(rast) gv
      FROM yourRasterTable
      WHERE rid = 1) foo;

You can also visualise individual pixel values by replacing ST_DumpAsPolygons() with ST_PixelAsPolygons().

Second ST_AsRaster(geom, template_raster) uses template_raster only for alignment, not for dimensions. No variant of ST_AsRaster() will let you determine the exact size of the resulting raster. It will always be of the same extent as the source geometry.

If you want your raster to be of a precise dimension you have to "burn" the result of ST_AsRaster() into en existing raster using ST_MapAlgebra(rast1, band1, rast2, band2, expression):

WITH raster AS (
  SELECT ST_SetSRID(ST_AddBand(ST_MakeEmptyRaster(200, 200, -1409340, 2138359, 4), '8BUI'::text), 92001) rast
), polyasraster AS (
  SELECT ST_AsRaster(geom, rast) polyrast
  FROM raster, (SELECT geom 
                FROM yourPolygonTable
                WHERE poly_id = 1234) foo
)
SELECT ST_MapAlgebra(rast, 1, polyrast, 1, '[rast1] + [rast2]') rast
FROM raster, polyasraster;

You can visualize the result in OpenJump with by adding the ST_PixelAsPolygons() part:

WITH raster AS (
  SELECT ST_SetSRID(ST_AddBand(ST_MakeEmptyRaster(200, 200, -1409340, 2138359, 4), '8BUI'::text), 92001) rast
), polyasraster AS (
  SELECT ST_AsRaster(geom, rast) polyrast
  FROM raster, (SELECT geom 
                FROM yourPolygonTable
                WHERE poly_id = 1234) foo
), finalraster AS ( 
SELECT ST_MapAlgebra(rast, 1, polyrast, 1, '[rast1] + [rast2]') rast
FROM raster, polyasraster
)
SELECT (gv).geom, (gv).val
FROM (SELECT ST_PixelAsPolygons(rast) gv
      FROM finalraster) foo;

Your Answer

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

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