2

I am running a simple zonal statistics query in PostGIS. My code is below

DROP TABLE IF EXISTS output_aa;
CREATE TABLE output_aa AS
    SELECT id, geom, (ST_SummaryStats(ST_Union(ST_Clip(rast,geom)),true)).mean
    FROM raster_aa, polys_aa
    GROUP BY id, geom; 

That works, although it does create a new table, rather than appending to the input polygon table.

But, if I have say 100 rasters in my database and 100 matching polygons tables, how would I best go about looping through all of them?

The raster names increment through a predefined sequence of letter pairs (e.g. AA, AB, AC, BA, BB etc) rather than numbers. The output can be either 100 output polygons tables or a single merged polygons table with all the output polygons in it. I'm still pretty new to PostGIS.


In response to the comments made, I have made an initial stab at the code and got part of the way there.

do
$$
declare
    rst record;
    out_tn text;
begin
    for rst in (select *
           from pg_catalog.pg_tables
           where tablename LIKE 'my_raster_%' ) loop 
        out_tn := 'results_' || rst.tablename;
        execute format ('drop table if exists %I ',out_tn);
        execute format ('create table %I as '
            'select id, geom, (ST_SummaryStats(ST_Union(ST_Clip(rast,geom)),true)).mean '
            'from %I, polys '
            'group by id, geom ', out_tn, rst.tablename);
    end loop;
end;
$$

This allows me to loop through all the rasters in the database which meet the requirement but only allows me to compare them to a single polygon table.

How do I extend this so that I can match my polygons tables with my rasters?

i.e. so that polys_aa queries my_raster_aa and polys_ab queries my_raster_ab etc.

I have some huge datasets to run, so I don't think it's sensible to try and put all the polygons in a single table.

Am I coming at this from the right perspective?

4
  • You'll need to create a function and use a loop. There's lots of examples in the documentation. You can search the pg_catalog to get table names select * from pg_catalog.pg_tables where tablename similar to ('raster_name_[A-Z]{2}'). I'm not sure if that's helpful or not. Or, to generate a series of suffixes, you could use with list as ( select chr(generate_series(65,90)) letter ) select t1.letter||t2.letter suffix from list t1, list t2 where t1.letter = t2.letter;
    – jbalk
    Nov 23 at 23:35
  • I'd start by learning how to create a function and look at examples of loops. You could generate the series of suffixes, then loop through them with a for loop. postgresqltutorial.com/plpgsql-for-loop
    – jbalk
    Nov 23 at 23:40
  • Thank you both, I've edited the original post as I've got part of the answer now, but not all of it.
    – Ddee
    Nov 24 at 16:19
  • When you say 'matching polygon tables', do you mean that if you have my_raster_AA, there is a corresponding my_polygons_AA?
    – jbalk
    Nov 29 at 18:13
0

After rereading your question, I realized that your suffix sequence isn't fixed (AA,BB,CC,DD,etc.). It will be easier to search for your raster prefix in the table catalog rather than trying to replicate the suffixes and tablenames. This code searches the catalog for your raster prefix ('my_raster_'), loops through the rasters finding the polygon table with the matching suffix and running stats on the raster:

create or replace function loop_stats()
returns void as $$
declare
    raster_prefix text := 'my_raster_';
    poly_prefix text := 'my_polys_';
    rst record;
    raster_suffix text;
    in_r text;
    in_p text;
    out_tn text;
begin 
    for rst in (select * from pg_catalog.pg_tables where tablename similar to raster_prefix || '[A-Z]{2}') loop
        in_r := rst.tablename;
        select replace(rst.tablename,raster_prefix,'') into raster_suffix;
        in_p := poly_prefix || raster_suffix;
        out_tn := 'results_' || in_r;
        execute 'drop table if exists "' || out_tn || '"';
        execute 'create table "' || out_tn || '" as select id, geom, (st_summarystats(st_union(st_clip(rast,geom)),true)).mean from "' || in_r || '" as r, "' || in_p || '" as p group by id, geom';
        raise notice '% processed using %......Results in %...', in_r, in_p, out_tn;
    end loop;
end;
$$ language plpgsql;

select loop_stats();

Previous Answer - Only works with fixed prefix sequence (AA,BB,CC,DD,etc.):

This function creates a temporary table. It creates the raster/polygon table name pairs by adding the raster and poly prefixes to the generated series of letter suffixes (AA, BB, etc.) and stores them in the temp table. It then loops through the name pairs, running stats on each pair and creating a results table for each raster. Adjust the start_char and end_char as needed, as well as the raster and poly prefixes.

create or replace function loop_stats()
returns void as $$
declare
    --ascii characters (A-Z = 65-90, a-z = 97-122)
    --A
    start_char int := 65;
    --Z
    end_char int := 90;
    raster_prefix text := 'my_raster_';
    poly_prefix text := 'my_polys_';
    rst record;
    in_r text;
    in_p text;
    out_tn text;
begin 
    drop table if exists temp_raster_polys;
    create temp table temp_raster_polys (raster_name text, poly_name text);
    execute 'insert into temp_raster_polys (raster_name, poly_name) with list as ( select chr(generate_series(' || start_char || ',' || end_char ||')) letter ) select ''' || raster_prefix || '''||t1.letter||t2.letter, ''' || poly_prefix || '''||t1.letter||t2.letter from list t1, list t2 where t1.letter = t2.letter';
    for rst in (select * from temp_raster_polys) loop
        in_r := rst.raster_name;
        in_p := rst.poly_name;
        out_tn := 'results_' || rst.raster_name;
        execute 'drop table if exists "' || out_tn || '"';
        execute 'create table "' || out_tn || '" as select id, geom, (st_summarystats(st_union(st_clip(rast,geom)),true)).mean from "' || in_r || '" as r, "' || in_p || '" as p group by id, geom';
        raise notice '% processed...Results in %...', in_r, out_tn;
    end loop;
end;
$$ language plpgsql;

select loop_stats();

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.