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) Commented Jul 30, 2019 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. Commented Jul 30, 2019 at 14:49

1 Answer 1

4

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

5
  • Brilliant, that's exactly what I needed. Seems like it's also a lot faster than my approach.
    – karpfen
    Commented Jul 31, 2019 at 8:57
  • I'm glad it`s still working. Which PostGIS version are you using? Commented Jul 31, 2019 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
    Commented Aug 1, 2019 at 7:53
  • I am struggling to run 'ST_ExtractToRaster' and getting the ERROR: function st_extracttoraster(raster, unknown, unknown) does not exist. Is there a way to setup this function? I am on postGIS 3.1
    – Shana
    Commented Jul 14, 2022 at 16:00
  • You have to install the PostGIS Addons script to define the function. Commented Jul 26, 2022 at 20:35

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.