3

I have a PostGIS table with Point geometries in it. I want to create a raster object where each cell value is the number of points within a given raster cell. Cells with zero points in it can be 0, Null or anything else that lets me differentiate it from a plausible count output.

What I did so far is to create a new raster (adapted from this tutorial):

CREATE TABLE public.sampleraster(
  rid SERIAL primary key, rast raster
);

INSERT INTO public.sampleraster(rast)
SELECT ST_AddBand(
  -- Make empty raster
  ST_MakeEmptyRaster(
      5000, -- Raster width x (in pixels)
      3000, -- Raster width y (in pixels)
      -75, -- Upper left X coordinate
      45, -- Upper left Y coordinate
      0.008333333333333, -- X Cell Resolution (in degrees)
      0.008333333333222, -- Y Cell Resolution (in degrees)
      0, -- X skew
      0, -- Y skew
      4326 -- SRID (WGS 84)
    ),
    ARRAY [
      ROW(
        1,  -- Band index: sets this as the first band
        '32BUI'::text,  -- Pixel Type (string rep of ST_BandPixelType types)
        NULL,  -- Initialized pixel value
        NULL   -- Nodata Value
      )
    ]::addbandarg[]
);

The Problem is now to fill the raster with count values. I'm counting the Points per cell like this

SELECT
round (st_x (geom) / 0.008333333333333) * 0.008333333333333 lon,
round (st_y (geom) / 0.008333333333222) * 0.008333333333222 lat,
count(*) ct
FROM sampledata
WHERE
st_x (geom) BETWEEN -75 AND -33.33333
AND 
st_y (geom) BETWEEN 10 AND 35

GROUP BY round (st_x (geom) / 0.008333333333333), round (st_y (geom) / 0.008333333333222)
ORDER BY round (st_x (geom) / 0.008333333333333), round (st_y (geom) / 0.008333333333222)

Now I'm struggling to get the values into the raster. I think I need to get the values in a 2D array (going from the tutorial), but I don't know how to go about this. Any ideas here?

If somebody is interested in recreating the issue, you can use this query to generate a bunch of random points in the study area:

CREATE TABLE sampledata AS
WITH x AS
    (select generate_series (1,1000),
     st_setSRID(st_makepoint(
        (random()*50)::float - 120,
        (random()*25)::float + 30),
     4326) geom,
     random()*100::int num_a)
SELECT geom FROM x;
2
  • 1
    if i look in postgis documentation, what you want to achieve seems to be possible with St_SetValues and a count of points at the given coordinate. so try to make a table with coordX, coordY, Value and then use St_setValues. Then, syntax would be St_SetValues(raster, nb_band (1), coordX, coordY, Value) – Maximilien jaffrès Jul 30 '19 at 13:56
  • Exactly as suggested by Maximilien. As a general rule, don't use zero in place of null, it is bad practice, and can lead to some nasty to find statistical bugs. – John Powell Jul 30 '19 at 14:49
3

Install the PostGIS Addons and have a look at ST_ExtractToRaster() with the COUNT_OF_POINTS method. Should go like this after you created a reference raster:

CREATE TABLE count_of_point_coverage AS
SELECT ST_ExtractToRaster(
         ST_AddBand(ST_MakeEmptyRaster(rast), '32BF'::text, -9999, -9999), 
         'public', 
         'sampledata', 
         'geom', 
          NULL, 
         'COUNT_OF_POINTS') rast
FROM sampleraster;

You need to have an index on the geometry column of your point table.

More details in the "Method 2" section of this post: http://geospatialelucubrations.blogspot.com/2014/05/a-guide-to-rasterization-of-vector.html

3
  • Brilliant, that's exactly what I needed. Seems like it's also a lot faster than my approach. – karpfen Jul 31 '19 at 8:57
  • I'm glad it`s still working. Which PostGIS version are you using? – Pierre Racine Jul 31 '19 at 16:20
  • I used it in two settings, both worked without issues: (1) Manjaro Linux, Postgres 11.4, PostGIS 2.5.2 r17328 (2) Windows Server 2012, Postgres 11.1, PostGIS 2.5.1 r17027 – karpfen Aug 1 '19 at 7:53

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.