8

I have a table of points and I'm trying to visualize their distributions. Each point has a POINT geometry (WGS84) and some numeric attributes: a, b, c.

| id | geom       | a  | b | c |
|-----------------|----|---|---|
| 1  | POINT(...  | 12 | 9 | 2 |
| etc ...                      |

The first thing I tried was to project to an equal area projection (EPSG: 3175), use ST_SnapToGrid and then group on the resulting point geometry.

SELECT 
  avg(a), avg(b), avg(c), 
  ST_Transform(ST_SnapToGrid(ST_Transform(geom, 3175), 50), 4326) AS bin
FROM points
GROUP BY bin

This groups the points into 50x50 meter bins with their aggregate values which I can import into qgis to play with. I optimized this by adding another column to the points table called bin_id which was a foreign key to a bins table:

| id | geom      |
|----|-----------|
| 1  | POINT(... |

Then I could write queries like this:

SELECT avg_a, avg_b, avg_c, geom
FROM (
  SELECT 
    avg(a) AS avg_a,
    avg(b) AS avg_b,
    avg(c) AS avg_c,
    bin_id
  FROM points
  GROUP BY bin_id
) AS tmp
JOIN bins ON tmp.cell_id = bins.id

This all works fine but i want to use Rasters instead. All the examples I've seen show how to import existing rasters or shapefiles into PostGIS. But I need to know how to generate one from data already in PostGIS.

Ideally I want a raster with 3 bands: avg_a, avg_b & avg_c

Can someone point me in the right direction?

EDIT 1: I'm not opposed to doing it outside of the database.
EDIT 2: Seems like I can use the GDAL API to do it.

1

1 Answer 1

4

We wrote a couple of articles about this you might find useful: http://www.bostongis.com/blog/index.php?/archives/211-Waiting-for-PostGIS-2.1-ST_SetValues-The-road-to-thematic-maps.html

http://www.bostongis.com/blog/index.php?/archives/222-Using-ST_ColorMap-for-Thematic-Maps.html

You'll want to check out in PostGIS manual -- ST_AsRaster, ST_DumpValues (and other pixel setter functions) http://postgis.net/docs/manual-2.1/RT_reference.html#Raster_Pixel_Accessors raster ST_Union -

3
  • Why is the scalex parameter to St_MakeEmptyRaster negative?
    – ilia choly
    Commented Mar 10, 2014 at 14:31
  • 1
    You mean the scale_y = -20000. Y direction in raster generally goes in opposite direction of the spatial reference space.
    – Regina Obe
    Commented Mar 10, 2014 at 18:29
  • Yes I meant scaley and that makes sense.
    – ilia choly
    Commented Mar 10, 2014 at 18:58

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.