Re: [postgis-users] Help with ST_Within query issue...

[Roxanne Reid-Bennett]

   ?modified | timestamp with time zone | |? |
>> 
>> gsw=# SELECT P.latitude, P.longitude, P.altitude, ST_AsText(P.point)
>> FROM Phenomena_v AS P
>> gsw-# WHERE
>> gsw-# P.acquisition_time >= '2023-01-08T20:14:43-0700'
>> gsw-# AND P.acquisition_time <= '2023-01-09T20:14:43-0700'
>> gsw-# AND ST_Within(ST_Transform(P.Point, 4326),
>> ST_GeomFromText('POLYGON((34.67010 -119.53618,
>> gsw'# 34.67010 -116.02055,
>> gsw'# 32.72798 -116.02055,
>> gsw'# 32.72798 -119.53618,
>> gsw'# 34.67010 -119.53618))',4326))
>> gsw-# = '1';
>> ERROR:? Unknown geometry type: 1025 - Point
>> CONTEXT:? parallel worker
>> gsw=#
>> 
>> Postgres 14 (There are points inside the polygon):
>> 
>> psql gsw
>> psql (14.6 (Ubuntu 14.6-1.pgdg22.04+1))
>> Type "help" for help.
>> 
>> gsw=# \dx
>> List of installed extensions
>>  ?? Name?? | Version | Schema?? | Description
>> 
>> --+-++-
>>  ?adminpack??? | 2.1 | pg_catalog | administrative functions
>> for PostgreSQL
>>  ?pg_freespacemap? | 1.2 | public | examine the free space
>> map (FSM)
>>  ?pgstattuple? | 1.5 | public | show tuple-level statistics
>>  ?plpgsql? | 1.0 | pg_catalog | PL/pgSQL procedural language
>>  ?postgis? | 3.2.1?? | public | PostGIS geometry,
>> geography, and raster spatial types and functions
>>  ?postgis_raster?? | 3.2.1?? | public | PostGIS raster types and
>> functions
>>  ?postgis_topology | 3.2.1?? | topology?? | PostGIS topology spatial
>> types and functions
>>  ?postgres_fdw | 1.1 | public | foreign-data wrapper for
>> remote PostgreSQL servers
>> (8 rows)
>> 
>> gsw=# \d Phenomena_v
>>  ??? View "public.phenomena_v"
>>  ??? Column??? | Type?? | Collation | Nullable |
>> Default
>> 
>> --+--+---+--+-
>>  ?id_pk??? | bigint?? |??
>> |? |
>>  ?meta_data_fk | bigint?? |??
>> |? |
>>  ?station_id?? | character varying(256)?? |??
>> |? |
>>  ?sensor_id??? | character varying(256)?? |??
>> |? |
>>  ?acquisition_time | timestamp with time zone |??
>> |? |
>>  ?acquisition_duration | bigint?? |??
>> |? |
>>  ?acquisition_period?? | bigint?? |??
>> |? |
>>  ?point??? | geometry(Point,4326) |??
>> |? |
>>  ?latitude | double precision |??
>> |? |
>>  ?longitude??? | double precision |??
>> |? |
>>  ?altitude | double precision |??
>> |? |
>>  ?dimension_id | character varying(256)?? |??
>> |? |
>>  ?quantity | character varying(255)?? |??
>> |? |
>>  ?units_id | character varying(256)?? |??
>> |? |
>>  ?created? | timestamp with time zone |??
>> |? |
>>  ?modified???????? | timestamp with time zone |??
>> |? |
>> 
>> gsw=# SELECT P.latitude, P.longitude, P.altitude, ST_AsText(P.point)
>> FROM Phenomena_v AS P
>> WHERE
>>  ??? P.acquisition_time >= '2023-01-08T20:14:43-0700'
>>  ??? AND P.acquisition_time <= '2023-01-09T20:14:43-0700'
>>  ??? AND ST_Within(ST_Transform(P.Point, 4326),
>> ST_GeomFromText('POLYGON((34.67010 -119.53618,
>> 34.67010 -116.02055,
>> 32.72798 -116.02055,
>> 32.72798 -119.53618,
>> 34.67010 -119.53618))',4326))
>> = '1';
>>  ?latitude | longitude | altitude | st_astext
>> --+---+--+---
>> (0 rows)
>> 
>> gsw=#
>> 
>> Thanks in advance.
>> 
>> --J
>> 
>> -- next part --
>&

Re: [postgis-users] Help with ST_Within query issue...

[Jeffrey Peacock]

  ?adminpack??? | 2.1 | pg_catalog | administrative functions
 for PostgreSQL
  ?pg_freespacemap? | 1.2 | public | examine the free space
 map (FSM)
  ?pgstattuple? | 1.5 | public | show tuple-level statistics
  ?plpgsql? | 1.0 | pg_catalog | PL/pgSQL procedural language
  ?postgis? | 3.2.1?? | public | PostGIS geometry,
 geography, and raster spatial types and functions
  ?postgis_raster?? | 3.2.1?? | public | PostGIS raster types and
 functions
  ?postgis_topology | 3.2.1?? | topology?? | PostGIS topology spatial
 types and functions
  ?postgres_fdw | 1.1 | public | foreign-data wrapper for
 remote PostgreSQL servers
 (8 rows)

 gsw=# \d Phenomena_v
  ??? View "public.phenomena_v"
  ??? Column??? | Type?? | Collation | Nullable |
 Default
 
--+--+---+--+-
  ?id_pk??? | bigint?? |??
 |? |
  ?meta_data_fk | bigint?? |??
 |? |
  ?station_id?? | character varying(256)?? |??
 |? |
  ?sensor_id??? | character varying(256)?? |??
 |? |
  ?acquisition_time | timestamp with time zone |??
 |? |
  ?acquisition_duration | bigint?? |??
 |? |
  ?acquisition_period?? | bigint?? |??
 |? |
  ?point??? | geometry(Point,4326) |??
 |? |
  ?latitude | double precision |??
 |? |
  ?longitude??? | double precision |??
 |? |
  ?altitude | double precision |??
 |? |
  ?dimension_id | character varying(256)?? |??
 |? |
  ?quantity | character varying(255)?? |??
 |? |
  ?units_id | character varying(256)?? |??
 |? |
  ?created? | timestamp with time zone |??
 |? |
  ?modified | timestamp with time zone |??
 |? |

 gsw=# SELECT P.latitude, P.longitude, P.altitude, ST_AsText(P.point)
 FROM Phenomena_v AS P
 WHERE
  ??? P.acquisition_time >= '2023-01-08T20:14:43-0700'
  ??? AND P.acquisition_time <= '2023-01-09T20:14:43-0700'
  ??? AND ST_Within(ST_Transform(P.Point, 4326),
 ST_GeomFromText('POLYGON((34.67010 -119.53618,
 34.67010 -116.02055,
 32.72798 -116.02055,
 32.72798 -119.53618,
 34.67010 -119.53618))',4326))
 = '1';
  ?latitude | longitude | altitude | st_astext
 --+---+--+---
 (0 rows)

 gsw=#

Thanks in advance.

--J

-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.osgeo.org/pipermail/postgis-users/attachments/20230111/40fbb315/attachment-0001.htm>

--

Message: 2
Date: Wed, 11 Jan 2023 08:35:59 -0800
From: Paul Ramsey 
To: PostGIS Users Discussion 
Subject: Re: [postgis-users] Help with ST_Within query issue...
Message-ID: 
Content-Type: text/plain;   charset=us-ascii

The error on pg10 is interesting and worrying, but at a first approximation 
your pg14 problem is that you have reversed the coordinate order in your 
polygon. The order should be longitude/latitude and yours is latitude/longitude.

Try

ST_GeomFromText('POLYGON((-119.53618 34.67010 ,
  -116.02055 34.67010 ,
  -116.02055 32.72798,
  -119.53618 32.72798 ,
  -119.53618 34.67010))', 4326)

P.


On Jan 11, 2023, at 8:31 AM, Jeffrey Peacock  wrote:

POLYGON((34.67010 -119.53618,



--

Message: 3
Date: Wed, 11 Jan 2023 11:40:52 -0500
From: "Regina Obe" 
To: "'PostGIS Users Discussion'" 
Subject: Re: [postgis-users] Help with ST_Within query issue...
Message-ID: <004801d925db$7866d150$693473f0$@pcorp.us>
Content-Type: text/plain; charset="utf-8"

Never seen that error before.  Does your view reference any foreign tables in 
another database.

My only guess is maybe some confusion with it reading the geometry table from 
another database.

  


Does a query like:

  


SELECT P.latitude, P.longitude, P.altitude, ST_AsText(P.point)

FROM Phenomena_v AS P

LIMIT 10;

  


If that doesn?t work what about

  


SELECT P.latitude, P.longitude, P.altitude

FROM Phenomena_v AS P

LIMIT 10;

  


If you can provide the definition of that view, that would be great.

  


Thanks,

Regina

  

  


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Jeffrey Peacock
Sent: Wednesday, January 11, 2023 11:32 AM
To:

Re: [postgis-users] Help with ST_Within query issue...

[Regina Obe]

Never seen that error before.  Does your view reference any foreign tables in 
another database.

My only guess is maybe some confusion with it reading the geometry table from 
another database.

 

Does a query like:

 

SELECT P.latitude, P.longitude, P.altitude, ST_AsText(P.point) 

FROM Phenomena_v AS P

LIMIT 10;

 

If that doesn’t work what about

 

SELECT P.latitude, P.longitude, P.altitude

FROM Phenomena_v AS P

LIMIT 10;

 

If you can provide the definition of that view, that would be great.

 

Thanks,

Regina

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Jeffrey Peacock
Sent: Wednesday, January 11, 2023 11:32 AM
To: postgis-users@lists.osgeo.org
Subject: [postgis-users] Help with ST_Within query issue...

 

Looking for some insight into why this is happening (2 examples below).

For Postgresql 10:

psql gsw
psql (14.6 (Ubuntu 14.6-1.pgdg20.04+1), server 10.23 (Ubuntu 
10.23-1.pgdg20.04+1))
Type "help" for help.

gsw=# \dx
 List of installed extensions
   Name   | Version |   Schema   | 
Description 
--+-++-
 plpgsql  | 1.0 | pg_catalog | PL/pgSQL procedural language
 postgis  | 3.2.1   | public | PostGIS geometry, geography, and 
raster spatial types and functions
 postgis_raster   | 3.2.1   | public | PostGIS raster types and functions
 postgis_topology | 3.2.1   | topology   | PostGIS topology spatial types and 
functions
 postgres_fdw | 1.0 | public | foreign-data wrapper for remote 
PostgreSQL servers
 tablefunc| 1.0 | public | functions that manipulate whole 
tables, including crosstab
 uuid-ossp| 1.0 | public | generate universally unique 
identifiers (UUIDs)
(7 rows)

gsw=# \d Phenomena_v
View "public.phenomena_v"
Column|   Type   | Collation | Nullable | 
Default 
--+--+---+--+-
 id_pk| bigint   |   |  | 
 meta_data_fk | bigint   |   |  | 
 station_id   | character varying(256)   |   |  | 
 sensor_id| character varying(256)   |   |  | 
 acquisition_time | timestamp with time zone |   |  | 
 acquisition_duration | bigint   |   |  | 
 acquisition_period   | bigint   |   |  | 
 point| geometry(Point,4326) |   |  | 
 latitude | double precision |   |  | 
 longitude| double precision |   |  | 
 altitude | double precision |   |  | 
 dimension_id | character varying(256)   |   |  | 
 quantity | character varying(255)   |   |  | 
 units_id | character varying(256)   |   |  | 
 created  | timestamp with time zone |   |  | 
 modified | timestamp with time zone |   |  | 

gsw=# SELECT P.latitude, P.longitude, P.altitude, ST_AsText(P.point) FROM 
Phenomena_v AS P
gsw-# WHERE
gsw-# P.acquisition_time >= '2023-01-08T20:14:43-0700'
gsw-# AND P.acquisition_time <= '2023-01-09T20:14:43-0700'
gsw-# AND ST_Within(ST_Transform(P.Point, 4326), 
ST_GeomFromText('POLYGON((34.67010 -119.53618,
gsw'#   
   34.67010 -116.02055,
gsw'#   
   32.72798 -116.02055,
gsw'#   
   32.72798 -119.53618,
gsw'#   
   34.67010 -119.53618))',4326))
gsw-#   
   = '1';
ERROR:  Unknown geometry type: 1025 - Point
CONTEXT:  parallel worker
gsw=# 

Postgres 14 (There are points inside the polygon):

psql gsw
psql (14.6 (Ubuntu 14.6-1.pgdg22.04+1))
Type "help" for help.

gsw=# \dx
 List of installed extensions
   Name   | Version |   Schema   | 
Description 
--+-++-
 adminpack| 2.1 | pg_catalog | administrative functions for 
PostgreSQL
 pg_freespacemap  | 1.2 | public | examine the free space map (FSM)
 pgstattuple  | 1.5 | public | show tuple-level statistics
 plpgsql  

Re: [postgis-users] Help with ST_Within query issue...

[Paul Ramsey]

The error on pg10 is interesting and worrying, but at a first approximation 
your pg14 problem is that you have reversed the coordinate order in your 
polygon. The order should be longitude/latitude and yours is latitude/longitude.

Try

ST_GeomFromText('POLYGON((-119.53618 34.67010 ,
 -116.02055 34.67010 ,
 -116.02055 32.72798,
 -119.53618 32.72798 ,
 -119.53618 34.67010))', 4326)

P.

> On Jan 11, 2023, at 8:31 AM, Jeffrey Peacock  wrote:
> 
> POLYGON((34.67010 -119.53618,

___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] help with a query

[Birgit Laggner]

Hi Jonatan, I am glad I could help. Have success with the tweaking!
Regards, Birgit

Am 26.09.2016 um 14:24 schrieb Jonatan Malaver:
Hi Birgit, thank you so much!! your code put me on the right track. 
I'm still tweaking it a bit since your code assumes only 1 line is 
attached to the endpoint. However, the endpoint of one line could 
break into 2 new starpoints.

thanks,
jon

On Tue, Sep 20, 2016 at 5:38 AM Birgit Laggner 
> wrote:


Hi Jonatan,

based on the anonymous code block Leknín came up with, I tried to
adapt it to your problem:

DO $$
DECLARE c_line record; c_gid integer; c_geom geometry; line
record; gid integer; geom geometry; i integer; n integer;
BEGIN

EXECUTE 'SELECT count(*) FROM electric_line' INTO n;

--initialize start
EXECUTE 'SELECT gid, geom FROM electric_line WHERE gid = 3312'
INTO c_line;

c_gid := c_line.gid;
c_geom := c_line.geom;
i := 1;

-- loop through lines following the flow direction
LOOP
EXECUTE 'SELECT gid, geom FROM electric_line
 WHERE
ST_DWithin(ST_EndPoint($1,ST_StartPoint(geom),0.01) OR
ST_DWithin(ST_EndPoint($1),ST_EndPoint(geom),0.01)' INTO line
USING c_geom;

gid := line.gid;
geom := line.geom;

--compare end point parent line with start point child line
and reverse line if necessary
IF NOT ST_DWithin(ST_EndPoint(c_geom),
ST_StartPoint(geom),0.01) THEN
EXECUTE
'UPDATE electric_line SET geom = ST_Reverse(geom)
WHERE gid = '||line.gid;
END IF;

--take child line as parent line for next loop
c_gid := gid;
c_geom := geom;
i := i + 1;

EXIT WHEN i = n;

END LOOP;

END $$;


Basically, you would start with your starting point for the flow.
Then you search for the next matching line and check if the
direction is ok. Otherwise you reverse the line. Then you go into
your next search loop and so on. The loop should exit when all
lines have been through the loop. I couldn't test the code but
hope you get the idea at least.

Regards,

Birgit



Am 14.09.2016 um 17:13 schrieb Jonatan Malaver:

Hello again, I do not have parent line id. All I have is a
starting point from where the direction should reference.

On Wed, Sep 14, 2016 at 9:09 AM Leknín Řepánek
> wrote:

On Wed, Sep 14, 2016 at 12:09:23PM +, Jonatan Malaver wrote:
> the reason being is that I do a network analysis by running
the following
> function:
> WITH RECURSIVE flow(gid, geom) AS (
> SELECT e.gid, e.geom FROM electric_line e, transformers
t WHERE ST_Distance
> (t.geom,ST_StartPoint(e.geom)) <= 0.01 AND t.gid=$1
>   UNION ALL
> SELECT n.gid, n.geom
> FROM electric_line n, flow f
> WHERE
ST_Distance(ST_EndPoint(f.geom),ST_StartPoint(n.geom)) <= 0.01
>   )
> The problem I have is that some of the lines direction are
in reversed. I'm
> trying to correct them with referenced to the first line.
Otherwise, I will end
> up changing hundreds of lines manually.
Manually? No. There are many of possible ways how do this by
query. For
example if you have in every line id of "parent line" you can use
anonymous block of code by something like this.
DO $$
DECLARE line record;
BEGIN
FOR line in SELECT lines from lines ORDER BY id LOOP
IF NOT ST_Equal(ST_startpoint(line.geom) , (SELECT
ST_EndPoint(geom) FROM lines WHERE id =
line.parent_line_id))
THEN
UPDATE lines SET geom = ST_Reverse(geom)
WHERE id =
line.id ;

END LOOP;

END

$$;


>
> On Tue, Sep 13, 2016 at 11:12 AM James Keener
> wrote:
>
> Depends on what you mean by direction. If you want to
grab the start and
> end points (st_startpoint and st_endpoint) and check
their x and y (st_x
> and st_y) for some condition (both less at the end?)
Then update the record
> with the value of st_reverse.
>
> I guess my other question is why it matters.
>
> Jim
>
> On September 13, 2016 8:31:07 AM EDT, Jonatan Malaver <
> jon.mala...@shrewsburyma.gov
> wrote:
>
> Hello,
>
>I'm trying to come up with a query that would
check the direction of
> a line. If the end point is 

Re: [postgis-users] help with a query

[Jonatan Malaver]

Hi Birgit, thank you so much!! your code put me on the right track. I'm
still tweaking it a bit since your code assumes only 1 line is attached to
the endpoint. However, the endpoint of one line could break into 2 new
starpoints.
thanks,
jon

On Tue, Sep 20, 2016 at 5:38 AM Birgit Laggner 
wrote:

> Hi Jonatan,
>
> based on the anonymous code block Leknín came up with, I tried to adapt it
> to your problem:
>
> DO $$
> DECLARE c_line record; c_gid integer; c_geom geometry; line record; gid
> integer; geom geometry; i integer; n integer;
> BEGIN
>
> EXECUTE 'SELECT count(*) FROM electric_line' INTO n;
>
> --initialize start
> EXECUTE 'SELECT gid, geom FROM electric_line WHERE gid = 3312' INTO c_line;
>
> c_gid := c_line.gid;
> c_geom := c_line.geom;
> i := 1;
>
> -- loop through lines following the flow direction
> LOOP
> EXECUTE 'SELECT gid, geom FROM electric_line
>  WHERE ST_DWithin(ST_EndPoint($1,ST_StartPoint(geom),0.01) OR
>ST_DWithin(ST_EndPoint($1),ST_EndPoint(geom),0.01)'
> INTO line USING c_geom;
>
> gid := line.gid;
> geom := line.geom;
>
> --compare end point parent line with start point child line and
> reverse line if necessary
> IF NOT ST_DWithin(ST_EndPoint(c_geom), ST_StartPoint(geom),0.01)
> THEN
> EXECUTE
> 'UPDATE electric_line SET geom = ST_Reverse(geom) WHERE
> gid = '||line.gid;
> END IF;
>
> --take child line as parent line for next loop
> c_gid := gid;
> c_geom := geom;
> i := i + 1;
>
> EXIT WHEN i = n;
>
> END LOOP;
>
> END $$;
>
>
> Basically, you would start with your starting point for the flow. Then you
> search for the next matching line and check if the direction is ok.
> Otherwise you reverse the line. Then you go into your next search loop and
> so on. The loop should exit when all lines have been through the loop. I
> couldn't test the code but hope you get the idea at least.
>
> Regards,
>
> Birgit
>
>
>
>  Am 14.09.2016 um 17:13 schrieb Jonatan Malaver:
>
> Hello again, I do not have parent line id. All I have is a starting point
> from where the direction should reference.
>
> On Wed, Sep 14, 2016 at 9:09 AM Leknín Řepánek 
> wrote:
>
>> On Wed, Sep 14, 2016 at 12:09:23PM +, Jonatan Malaver wrote:
>> > the reason being is that I do a network analysis by running the
>> following
>> > function:
>> > WITH RECURSIVE flow(gid, geom) AS (
>> > SELECT e.gid, e.geom FROM electric_line e, transformers t WHERE
>> ST_Distance
>> > (t.geom,ST_StartPoint(e.geom)) <= 0.01 AND t.gid=$1
>> >   UNION ALL
>> > SELECT n.gid, n.geom
>> > FROM electric_line n, flow f
>> > WHERE ST_Distance(ST_EndPoint(f.geom),ST_StartPoint(n.geom)) <= 0.01
>> >   )
>> > The problem I have is that some of the lines direction are in reversed.
>> I'm
>> > trying to correct them with referenced to the first line. Otherwise, I
>> will end
>> > up changing hundreds of lines manually.
>> Manually? No. There are many of possible ways how do this by query. For
>> example if you have in every line id of "parent line" you can use
>> anonymous block of code by something like this.
>> DO $$
>> DECLARE line record;
>> BEGIN
>> FOR line in SELECT lines from lines ORDER BY id LOOP
>> IF NOT ST_Equal(ST_startpoint(line.geom) , (SELECT
>> ST_EndPoint(geom) FROM lines WHERE id = line.parent_line_id))
>> THEN
>> UPDATE lines SET geom = ST_Reverse(geom) WHERE id =
>> line.id;
>>
>> END LOOP;
>>
>> END
>>
>> $$;
>>
>>
>> >
>> > On Tue, Sep 13, 2016 at 11:12 AM James Keener 
>> wrote:
>> >
>> > Depends on what you mean by direction. If you want to grab the
>> start and
>> > end points (st_startpoint and st_endpoint) and check their x and y
>> (st_x
>> > and st_y) for some condition (both less at the end?) Then update
>> the record
>> > with the value of st_reverse.
>> >
>> > I guess my other question is why it matters.
>> >
>> > Jim
>> >
>> > On September 13, 2016 8:31:07 AM EDT, Jonatan Malaver <
>> > jon.mala...@shrewsburyma.gov> wrote:
>> >
>> > Hello,
>> >
>> >I'm trying to come up with a query that would check the
>> direction of
>> > a line. If the end point is not the start point of the next
>> line to
>> > update the line by reversing that line. Can anyone give me
>> pointers on
>> > how to do it?
>> >
>> > Thanks,
>> > Jon
>> >
>> >
>> > postgis-users mailing list
>> > postgis-users@lists.osgeo.org
>> > http://lists.osgeo.org/mailman/listinfo/postgis-users
>> >
>> >
>> > --
>> > Sent from my Android device with K-9 Mail. Please excuse my brevity.
>> >
>> > --
>> >
>> > Thanks,
>> >
>> >
>> > Jonatan Malaver
>> >
>> > Assistant Engineer of Electrical and Cable Operations
>> >
>> > Shrewsbury Electric & Cable Operations
>> >
>> > 100 

Re: [postgis-users] help with a query

[Birgit Laggner]

Hi Jonatan,

based on the anonymous code block Leknín came up with, I tried to adapt 
it to your problem:


DO $$
DECLARE c_line record; c_gid integer; c_geom geometry; line record; gid 
integer; geom geometry; i integer; n integer;

BEGIN

EXECUTE 'SELECT count(*) FROM electric_line' INTO n;

--initialize start
EXECUTE 'SELECT gid, geom FROM electric_line WHERE gid = 3312' INTO c_line;

c_gid := c_line.gid;
c_geom := c_line.geom;
i := 1;

-- loop through lines following the flow direction
LOOP
EXECUTE 'SELECT gid, geom FROM electric_line
 WHERE ST_DWithin(ST_EndPoint($1,ST_StartPoint(geom),0.01) OR
ST_DWithin(ST_EndPoint($1),ST_EndPoint(geom),0.01)' INTO line USING c_geom;

gid := line.gid;
geom := line.geom;

--compare end point parent line with start point child line and 
reverse line if necessary
IF NOT ST_DWithin(ST_EndPoint(c_geom), 
ST_StartPoint(geom),0.01) THEN

EXECUTE
'UPDATE electric_line SET geom = ST_Reverse(geom) WHERE 
gid = '||line.gid;

END IF;

--take child line as parent line for next loop
c_gid := gid;
c_geom := geom;
i := i + 1;

EXIT WHEN i = n;

END LOOP;

END $$;


Basically, you would start with your starting point for the flow. Then 
you search for the next matching line and check if the direction is ok. 
Otherwise you reverse the line. Then you go into your next search loop 
and so on. The loop should exit when all lines have been through the 
loop. I couldn't test the code but hope you get the idea at least.


Regards,

Birgit


Am 14.09.2016 um 17:13 schrieb Jonatan Malaver:
Hello again, I do not have parent line id. All I have is a starting 
point from where the direction should reference.


On Wed, Sep 14, 2016 at 9:09 AM Leknín Řepánek 
> wrote:


On Wed, Sep 14, 2016 at 12:09:23PM +, Jonatan Malaver wrote:
> the reason being is that I do a network analysis by running the
following
> function:
> WITH RECURSIVE flow(gid, geom) AS (
> SELECT e.gid, e.geom FROM electric_line e, transformers t
WHERE ST_Distance
> (t.geom,ST_StartPoint(e.geom)) <= 0.01 AND t.gid=$1
>   UNION ALL
> SELECT n.gid, n.geom
> FROM electric_line n, flow f
> WHERE ST_Distance(ST_EndPoint(f.geom),ST_StartPoint(n.geom))
<= 0.01
>   )
> The problem I have is that some of the lines direction are in
reversed. I'm
> trying to correct them with referenced to the first line.
Otherwise, I will end
> up changing hundreds of lines manually.
Manually? No. There are many of possible ways how do this by
query. For
example if you have in every line id of "parent line" you can use
anonymous block of code by something like this.
DO $$
DECLARE line record;
BEGIN
FOR line in SELECT lines from lines ORDER BY id LOOP
IF NOT ST_Equal(ST_startpoint(line.geom) , (SELECT
ST_EndPoint(geom) FROM lines WHERE id = line.parent_line_id))
THEN
UPDATE lines SET geom = ST_Reverse(geom) WHERE id =
line.id ;

END LOOP;

END

$$;


>
> On Tue, Sep 13, 2016 at 11:12 AM James Keener > wrote:
>
> Depends on what you mean by direction. If you want to grab
the start and
> end points (st_startpoint and st_endpoint) and check their x
and y (st_x
> and st_y) for some condition (both less at the end?) Then
update the record
> with the value of st_reverse.
>
> I guess my other question is why it matters.
>
> Jim
>
> On September 13, 2016 8:31:07 AM EDT, Jonatan Malaver <
> jon.mala...@shrewsburyma.gov
> wrote:
>
> Hello,
>
>I'm trying to come up with a query that would check
the direction of
> a line. If the end point is not the start point of the
next line to
> update the line by reversing that line. Can anyone give
me pointers on
> how to do it?
>
> Thanks,
> Jon
>
>
> postgis-users mailing list
> postgis-users@lists.osgeo.org 
> http://lists.osgeo.org/mailman/listinfo/postgis-users
>
>
> --
> Sent from my Android device with K-9 Mail. Please excuse my
brevity.
>
> --
>
> Thanks,
>
>
> Jonatan Malaver
>
> Assistant Engineer of Electrical and Cable Operations
>
> Shrewsbury Electric & Cable Operations
>
> 100 Maple Avenue
>
> Shrewsbury, MA 01545
>
> Office: (508) 841-8610
>
> Fax: (508) 842-9267
>
> jon.mala...@shrewsburyma.gov 
>

> 

Re: [postgis-users] help with a query

[Jonatan Malaver]

Hello again, I do not have parent line id. All I have is a starting point
from where the direction should reference.

On Wed, Sep 14, 2016 at 9:09 AM Leknín Řepánek 
wrote:

> On Wed, Sep 14, 2016 at 12:09:23PM +, Jonatan Malaver wrote:
> > the reason being is that I do a network analysis by running the following
> > function:
> > WITH RECURSIVE flow(gid, geom) AS (
> > SELECT e.gid, e.geom FROM electric_line e, transformers t WHERE
> ST_Distance
> > (t.geom,ST_StartPoint(e.geom)) <= 0.01 AND t.gid=$1
> >   UNION ALL
> > SELECT n.gid, n.geom
> > FROM electric_line n, flow f
> > WHERE ST_Distance(ST_EndPoint(f.geom),ST_StartPoint(n.geom)) <= 0.01
> >   )
> > The problem I have is that some of the lines direction are in reversed.
> I'm
> > trying to correct them with referenced to the first line. Otherwise, I
> will end
> > up changing hundreds of lines manually.
> Manually? No. There are many of possible ways how do this by query. For
> example if you have in every line id of "parent line" you can use
> anonymous block of code by something like this.
> DO $$
> DECLARE line record;
> BEGIN
> FOR line in SELECT lines from lines ORDER BY id LOOP
> IF NOT ST_Equal(ST_startpoint(line.geom) , (SELECT
> ST_EndPoint(geom) FROM lines WHERE id = line.parent_line_id))
> THEN
> UPDATE lines SET geom = ST_Reverse(geom) WHERE id =
> line.id;
>
> END LOOP;
>
> END
>
> $$;
>
>
> >
> > On Tue, Sep 13, 2016 at 11:12 AM James Keener  wrote:
> >
> > Depends on what you mean by direction. If you want to grab the start
> and
> > end points (st_startpoint and st_endpoint) and check their x and y
> (st_x
> > and st_y) for some condition (both less at the end?) Then update the
> record
> > with the value of st_reverse.
> >
> > I guess my other question is why it matters.
> >
> > Jim
> >
> > On September 13, 2016 8:31:07 AM EDT, Jonatan Malaver <
> > jon.mala...@shrewsburyma.gov> wrote:
> >
> > Hello,
> >
> >I'm trying to come up with a query that would check the
> direction of
> > a line. If the end point is not the start point of the next line
> to
> > update the line by reversing that line. Can anyone give me
> pointers on
> > how to do it?
> >
> > Thanks,
> > Jon
> >
> >
> > postgis-users mailing list
> > postgis-users@lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/postgis-users
> >
> >
> > --
> > Sent from my Android device with K-9 Mail. Please excuse my brevity.
> >
> > --
> >
> > Thanks,
> >
> >
> > Jonatan Malaver
> >
> > Assistant Engineer of Electrical and Cable Operations
> >
> > Shrewsbury Electric & Cable Operations
> >
> > 100 Maple Avenue
> >
> > Shrewsbury, MA 01545
> >
> > Office: (508) 841-8610
> >
> > Fax: (508) 842-9267
> >
> > jon.mala...@shrewsburyma.gov
> >
>
> > ___
> > postgis-users mailing list
> > postgis-users@lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/postgis-users
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] help with a query

[Leknín Řepánek]

On Wed, Sep 14, 2016 at 12:09:23PM +, Jonatan Malaver wrote:
> the reason being is that I do a network analysis by running the following
> function:
> WITH RECURSIVE flow(gid, geom) AS (
>     SELECT e.gid, e.geom FROM electric_line e, transformers t WHERE 
> ST_Distance
> (t.geom,ST_StartPoint(e.geom)) <= 0.01 AND t.gid=$1
>   UNION ALL
>     SELECT n.gid, n.geom
>     FROM electric_line n, flow f
>     WHERE ST_Distance(ST_EndPoint(f.geom),ST_StartPoint(n.geom)) <= 0.01
>   )
> The problem I have is that some of the lines direction are in reversed. I'm
> trying to correct them with referenced to the first line. Otherwise, I will 
> end
> up changing hundreds of lines manually.
Manually? No. There are many of possible ways how do this by query. For
example if you have in every line id of "parent line" you can use
anonymous block of code by something like this.
DO $$
DECLARE line record;
BEGIN
FOR line in SELECT lines from lines ORDER BY id LOOP
IF NOT ST_Equal(ST_startpoint(line.geom) , (SELECT
ST_EndPoint(geom) FROM lines WHERE id = line.parent_line_id))
THEN
UPDATE lines SET geom = ST_Reverse(geom) WHERE id =
line.id;

END LOOP;

END

$$;


> 
> On Tue, Sep 13, 2016 at 11:12 AM James Keener  wrote:
> 
> Depends on what you mean by direction. If you want to grab the start and
> end points (st_startpoint and st_endpoint) and check their x and y (st_x
> and st_y) for some condition (both less at the end?) Then update the 
> record
> with the value of st_reverse.
> 
> I guess my other question is why it matters.
> 
> Jim
> 
> On September 13, 2016 8:31:07 AM EDT, Jonatan Malaver <
> jon.mala...@shrewsburyma.gov> wrote:
> 
> Hello,
> 
>    I'm trying to come up with a query that would check the direction 
> of
> a line. If the end point is not the start point of the next line to
> update the line by reversing that line. Can anyone give me pointers on
> how to do it?
> 
> Thanks,
> Jon
> 
> 
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/postgis-users
> 
> 
> --
> Sent from my Android device with K-9 Mail. Please excuse my brevity.
> 
> --
> 
> Thanks,
> 
> 
> Jonatan Malaver
> 
> Assistant Engineer of Electrical and Cable Operations
> 
> Shrewsbury Electric & Cable Operations
> 
> 100 Maple Avenue
> 
> Shrewsbury, MA 01545
> 
> Office: (508) 841-8610
> 
> Fax: (508) 842-9267
> 
> jon.mala...@shrewsburyma.gov
> 

> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] help with a query

[Jonatan Malaver]

the reason being is that I do a network analysis by running the following
function:
WITH RECURSIVE flow(gid, geom) AS (
SELECT e.gid, e.geom FROM electric_line e, transformers t WHERE
ST_Distance(t.geom,ST_StartPoint(e.geom)) <= 0.01 AND t.gid=$1
  UNION ALL
SELECT n.gid, n.geom
FROM electric_line n, flow f
WHERE ST_Distance(ST_EndPoint(f.geom),ST_StartPoint(n.geom)) <= 0.01
  )
The problem I have is that some of the lines direction are in reversed. I'm
trying to correct them with referenced to the first line. Otherwise, I will
end up changing hundreds of lines manually.

On Tue, Sep 13, 2016 at 11:12 AM James Keener  wrote:

> Depends on what you mean by direction. If you want to grab the start and
> end points (st_startpoint and st_endpoint) and check their x and y (st_x
> and st_y) for some condition (both less at the end?) Then update the record
> with the value of st_reverse.
>
> I guess my other question is why it matters.
>
> Jim
>
> On September 13, 2016 8:31:07 AM EDT, Jonatan Malaver <
> jon.mala...@shrewsburyma.gov> wrote:
>
>> Hello,
>>
>>I'm trying to come up with a query that would check the direction of a
>> line. If the end point is not the start point of the next line to update
>> the line by reversing that line. Can anyone give me pointers on how to do
>> it?
>>
>> Thanks,
>> Jon
>>
>>
>> postgis-users mailing list
>> postgis-users@lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/postgis-users
>>
>>
> --
> Sent from my Android device with K-9 Mail. Please excuse my brevity.
>
-- 

Thanks,


Jonatan Malaver

Assistant Engineer of Electrical and Cable Operations

Shrewsbury Electric & Cable Operations

100 Maple Avenue

Shrewsbury, MA 01545

Office: (508) 841-8610

Fax: (508) 842-9267

jon.mala...@shrewsburyma.gov
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] help with a query

[Jonatan Malaver]

This is the query I have so far:
WITH RECURSIVE flow(gid, geom) AS (
SELECT e.gid, e.geom FROM electric_line e WHERE e.gid = 3312
  UNION ALL
SELECT n.gid, n.geom
FROM electric_line n, flow f
WHERE ST_DWithin(ST_EndPoint(f.geom),ST_StartPoint(n.geom),0.01)
  )
UPDATE electric_line SET geom = ST_Reverse(electric_line.geom) FROM flow
WHERE electric_line.gid = flow.gid;

Is it possible to add an If statement that would only update if the
EndPoint is not StartPoint of the next the line?

On Tue, Sep 13, 2016 at 8:37 AM Leknín Řepánek 
wrote:

> ST_EndPoint, ST_StartPoint, ST_Reverse.
>
> On Tue, Sep 13, 2016 at 12:31:07PM +, Jonatan Malaver wrote:
> > Hello,
> >
> >I'm trying to come up with a query that would check the direction of
> a line.
> > If the end point is not the start point of the next line to update the
> line by
> > reversing that line. Can anyone give me pointers on how to do it?
> >
> > Thanks,
> > Jon
>
> > ___
> > postgis-users mailing list
> > postgis-users@lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/postgis-users
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] help with a query

[Leknín Řepánek]

ST_EndPoint, ST_StartPoint, ST_Reverse.

On Tue, Sep 13, 2016 at 12:31:07PM +, Jonatan Malaver wrote:
> Hello,
> 
>    I'm trying to come up with a query that would check the direction of a 
> line.
> If the end point is not the start point of the next line to update the line by
> reversing that line. Can anyone give me pointers on how to do it?
> 
> Thanks,
> Jon

> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Paragon Corporation]

Darrel,

 

Don't be quiet.  People rarely ask me raster questions, so it's a refreshing 
change :) to be able to sharpen my raster chops.

 

Thanks,

Regina

http://www.postgis.us

http://postgis.net

 

 

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Sunday, November 29, 2015 7:38 PM
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

 

Dear Regina,

 

Just to wrap this one up – you were right! 

 

I decided to make the network file in arcgis and then import that into my 
postgis db.  I then ran the same query but with the new polygon and the totals 
were exactly the same as the point file and the data I extracted in QGIS.

 

It is clear that I need to look at how I created the polygon in postgis 
(although whatever problems this caused they were ignored by QGIS).  Thankfully 
I can do what I need to do now (the polygon extraction takes around 25 minutes 
which is half the time the point file took and that will do !), and I am happy 
with the workflow. 

 

Many thanks again for all of your help and there is no need to reply. I promise 
to be quiet for a while now.

 

Best wishes

 

Darrel 

 

 

 

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 29 November, 2015 4:08 PM
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org 
<mailto:postgis-users@lists.osgeo.org> >
Subject: Re: [postgis-users] Help with SQL query?

 

Darrel,

 

I'm a bit confused about something.  I thought concentrated was what you were 
using to determine whether or not to consider a pixel in the deposition raster, 
so shouldn't the query be:

 

WITH  foo AS (

  SELECT  mymodel.deposition.rid,  mymodel.networkpoly.gid,  geom,  
ST_SummaryStats( ST_Clip(rast, geom) ) As st

FROM mymodel.deposition INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )

   WHERE filename='10_inci.tif'

)

SELECT sum( (st).sum)

FROM foo;

 

 

As your original query – was using deposition.

 

 

CREATE TABLE mymodel.networkdep AS

 

SELECT filename, gid, ST_Value(rast, geom) val

 

FROM mymodel.deposition, mymodel.network

 

WHERE ST_Intersects(rast, geom) 

 

ORDER BY gid, rid;

 

 

Gathering from the count it looks like all the pixels of those are treated as 
data when some should be treated as no data.  If it's not an error in the 
raster you are using,

 

it would really help to get output of the polygon and also the other stats I 
mentioned which (st).* outputs  So a query something like:

 

WITH  foo AS (

  SELECT  mymodel.concentrated.rid,  mymodel.networkpoly.gid,  geom,  
ST_SummaryStats( ST_Clip(rast, geom) ) As st

FROM mymodel.concentrated INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )

   WHERE filename='10_inci.tif'

)

SELECT rid, gid, geom, (st).*

FROM foo;

 

 

You might actually want to plot the geom on the map overlaid with the problem 
tiles (or original raster) you classified for that region.

 

My guess is something went wrong in the ST_Reclass and all the pixels in those 
tiles did not get set to 0 (and thus are treated as data doing a summary stats 
of the problem pre-reclassed tiles and looking at the min,max will give you a 
clue about that).

This second point may be moot if it's just you meant to put deposition instead 
of concentration, but good knowledge to keep in mind for future troubleshooting.

 

Hope that helps,

Regina

 

 

 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Dear Regina and Roxanne,

Many thanks for your suggestions.  I reran the query to show the breakdown of 
the sum elements as suggested.

Now I know (or at least think I know) from the points version that there should 
be around 15500 intersect,  Looking at the table of information there are some 
odd counts for two rids.

rid   countsum
1779   6553616.78682
1810   655363.576006

There are 186 rids in total and although there are duplicate rids there are no 
duplicate counts that I can see, so I think that is OK. The summation of the 
sum column gives the high value previously reported.

The counts on the two rids shown above make little sense and they do, of 
course, seem to account for the discrepancy in sums (albeit removing these 
would overshoot the alternate sum at 0.898929)

I’m afraid this does not make me any wiser, especially as if I load the polygon 
and raster into QGIS from the same postgis tables and do the extraction there ( 
using grid statistics for polygons and then save the attribute table to a csv), 
the summation  gives the lower sum and a cell count of 15k which looks right.

I must confess I have no idea why the sql query returns these unusual counts 
and sums. The actual raw data tables appear to be fine.

Darrel


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 29 November 2015 03:30
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,

Try not summing.  That might give you more information about  what is wrong.  
So instead do this:

WITH  foo AS (
  SELECT  mymodel.concentrated.rid,  mymodel.networkpoly.gid,  ST_SummaryStats( 
ST_Clip(rast, geom) ) As st
FROM mymodel.concentrated INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )
   WHERE filename='10_inci.tif'
)
SELECT rid, gid, (st).*
FROM foo;

That might give you multiple records per rid which might be okay since a single 
tile might be clipped by more than one polygon.  At anyrate, you'll get a lot 
more stats and be able to target the tile/poly combo that might be at issue.

So your query should have

rid, gid, count, sum, mean, stddev, min, max  (for each tile/geom combo – which 
should hopefully give you a clue what went wrong).

Hope that helps,
Regina



___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Paragon Corporation]

Darrel,

 

I'm a bit confused about something.  I thought concentrated was what you were 
using to determine whether or not to consider a pixel in the deposition raster, 
so shouldn't the query be:

 

WITH  foo AS (

  SELECT  mymodel.deposition.rid,  mymodel.networkpoly.gid,  geom,  
ST_SummaryStats( ST_Clip(rast, geom) ) As st

FROM mymodel.deposition INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )

   WHERE filename='10_inci.tif'

)

SELECT sum( (st).sum)

FROM foo;

 

 

As your original query – was using deposition.

 

 

CREATE TABLE mymodel.networkdep AS

 

SELECT filename, gid, ST_Value(rast, geom) val

 

FROM mymodel.deposition, mymodel.network

 

WHERE ST_Intersects(rast, geom) 

 

ORDER BY gid, rid;

 

 

Gathering from the count it looks like all the pixels of those are treated as 
data when some should be treated as no data.  If it's not an error in the 
raster you are using,

 

it would really help to get output of the polygon and also the other stats I 
mentioned which (st).* outputs  So a query something like:

 

WITH  foo AS (

  SELECT  mymodel.concentrated.rid,  mymodel.networkpoly.gid,  geom,  
ST_SummaryStats( ST_Clip(rast, geom) ) As st

FROM mymodel.concentrated INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )

   WHERE filename='10_inci.tif'

)

SELECT rid, gid, geom, (st).*

FROM foo;

 

 

You might actually want to plot the geom on the map overlaid with the problem 
tiles (or original raster) you classified for that region.

 

My guess is something went wrong in the ST_Reclass and all the pixels in those 
tiles did not get set to 0 (and thus are treated as data doing a summary stats 
of the problem pre-reclassed tiles and looking at the min,max will give you a 
clue about that).

This second point may be moot if it's just you meant to put deposition instead 
of concentration, but good knowledge to keep in mind for future troubleshooting.

 

Hope that helps,

Regina

 

 

 

 

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Sunday, November 29, 2015 10:43 AM
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

 

Dear Regina and Roxanne,

 

Many thanks for your suggestions.  I reran the query to show the breakdown of 
the sum elements as suggested. 

 

Now I know (or at least think I know) from the points version that there should 
be around 15500 intersect,  Looking at the table of information there are some 
odd counts for two rids.

 

rid   countsum

1779   6553616.78682

1810   655363.576006

 

There are 186 rids in total and although there are duplicate rids there are no 
duplicate counts that I can see, so I think that is OK. The summation of the 
sum column gives the high value previously reported.

 

The counts on the two rids shown above make little sense and they do, of 
course, seem to account for the discrepancy in sums (albeit removing these 
would overshoot the alternate sum at 0.898929)

 

I’m afraid this does not make me any wiser, especially as if I load the polygon 
and raster into QGIS from the same postgis tables and do the extraction there ( 
using grid statistics for polygons and then save the attribute table to a csv), 
the summation  gives the lower sum and a cell count of 15k which looks right.

 

I must confess I have no idea why the sql query returns these unusual counts 
and sums. The actual raw data tables appear to be fine.

 

Darrel

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 29 November 2015 03:30
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org 
<mailto:postgis-users@lists.osgeo.org> >
Subject: Re: [postgis-users] Help with SQL query?

 

Darrel,

 

Try not summing.  That might give you more information about  what is wrong.  
So instead do this:

 

WITH  foo AS (

  SELECT  mymodel.concentrated.rid,  mymodel.networkpoly.gid,  ST_SummaryStats( 
ST_Clip(rast, geom) ) As st

FROM mymodel.concentrated INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )

   WHERE filename='10_inci.tif'

)

SELECT rid, gid, (st).*

FROM foo;

 

That might give you multiple records per rid which might be okay since a single 
tile might be clipped by more than one polygon.  At anyrate, you'll get a lot 
more stats and be able to target the tile/poly combo that might be at issue.

 

So your query should have

 

rid, gid, count, sum, mean, stddev, min, max  (for each tile/geom combo – which 
should hopefully give you a clue what went wrong).

 

Hope that helps,

Regina

 

 

 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Dear Regina,

Apologies I pasted from the same operation on a different table.

This operation is being used to extract various attributes from the same 
location. The data are stored in a number of tables (e.g. table 
concentrated[has rasters with filenames such as 10_inci.tif], table diffuse 
[filenames are 10_ero.tif], deposition [filenames are 10_depo] ..etc). There 
are 8 such tables – each have 136 tifs tiles at 256x256).  I’m afraid the file 
naming system is  little of a dogs breakfast – sorry.

The polygon for the network is derived from a flow accumulation raster 
(10_FA.tif etc) – I did this separately and saved the result as networkpoly 
which is used in all cases to extract from the various rasters.

Anyhow the values I used are from the concentrated table i.e. if I do this 
using intersection with the points file I get  a sum of ~0.9, I get the same 
sum value if I do the extraction in QGIS using networkpoly. If I use the query 
in postgres I get ~21 for the summation.

The networkpolygon looks correct in QGIS so I am a little confused over the 
suggestion that ST_Reclass may be the problem – that seems to have worked fine, 
there is certainly nothing visable that would account for intersection with an 
additional ~130k cells in the underlying raster. Anyhow I will continue to dig, 
but at present it is difficult to see why these would work in QGIS but not in 
postgres/postgis

Apologies again for the confusion with the filenames – the operation remains 
the same.

Darrel

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 29 November, 2015 4:08 PM
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,

I'm a bit confused about something.  I thought concentrated was what you were 
using to determine whether or not to consider a pixel in the deposition raster, 
so shouldn't the query be:

WITH  foo AS (
  SELECT  mymodel.deposition.rid,  mymodel.networkpoly.gid,  geom,  
ST_SummaryStats( ST_Clip(rast, geom) ) As st
FROM mymodel.deposition INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )
   WHERE filename='10_inci.tif'
)
SELECT sum( (st).sum)
FROM foo;


As your original query – was using deposition.


CREATE TABLE mymodel.networkdep AS

SELECT filename, gid, ST_Value(rast, geom) val

FROM mymodel.deposition, mymodel.network

WHERE ST_Intersects(rast, geom)

ORDER BY gid, rid;


Gathering from the count it looks like all the pixels of those are treated as 
data when some should be treated as no data.  If it's not an error in the 
raster you are using,

it would really help to get output of the polygon and also the other stats I 
mentioned which (st).* outputs  So a query something like:

WITH  foo AS (
  SELECT  mymodel.concentrated.rid,  mymodel.networkpoly.gid,  geom,  
ST_SummaryStats( ST_Clip(rast, geom) ) As st
FROM mymodel.concentrated INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )
   WHERE filename='10_inci.tif'
)
SELECT rid, gid, geom, (st).*
FROM foo;


You might actually want to plot the geom on the map overlaid with the problem 
tiles (or original raster) you classified for that region.

My guess is something went wrong in the ST_Reclass and all the pixels in those 
tiles did not get set to 0 (and thus are treated as data doing a summary stats 
of the problem pre-reclassed tiles and looking at the min,max will give you a 
clue about that).
This second point may be moot if it's just you meant to put deposition instead 
of concentration, but good knowledge to keep in mind for future troubleshooting.

Hope that helps,
Regina






From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Sunday, November 29, 2015 10:43 AM
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Dear Regina and Roxanne,

Many thanks for your suggestions.  I reran the query to show the breakdown of 
the sum elements as suggested.

Now I know (or at least think I know) from the points version that there should 
be around 15500 intersect,  Looking at the table of information there are some 
odd counts for two rids.

rid   countsum
1779   6553616.78682
1810   655363.576006

There are 186 rids in total and although there are duplicate rids there are no 
duplicate counts that I can see, so I think that is OK. The summation of the 
sum column gives the high value previously reported.

The counts on the two rids shown above make little sense and they do, of 
course, seem to account for the discrepancy in sums (albeit removing these 
would overshoot the alternate sum at 0.898929)

I’m afraid this does not make me any wiser, especially as if I load the polygon 
and raster into QGIS from the same postgis tables and do the extraction t

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Dear Regina,

Just to wrap this one up – you were right!

I decided to make the network file in arcgis and then import that into my 
postgis db.  I then ran the same query but with the new polygon and the totals 
were exactly the same as the point file and the data I extracted in QGIS.

It is clear that I need to look at how I created the polygon in postgis 
(although whatever problems this caused they were ignored by QGIS).  Thankfully 
I can do what I need to do now (the polygon extraction takes around 25 minutes 
which is half the time the point file took and that will do !), and I am happy 
with the workflow.

Many thanks again for all of your help and there is no need to reply. I promise 
to be quiet for a while now.

Best wishes

Darrel





From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 29 November, 2015 4:08 PM
To: 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,

I'm a bit confused about something.  I thought concentrated was what you were 
using to determine whether or not to consider a pixel in the deposition raster, 
so shouldn't the query be:

WITH  foo AS (
  SELECT  mymodel.deposition.rid,  mymodel.networkpoly.gid,  geom,  
ST_SummaryStats( ST_Clip(rast, geom) ) As st
FROM mymodel.deposition INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )
   WHERE filename='10_inci.tif'
)
SELECT sum( (st).sum)
FROM foo;


As your original query – was using deposition.


CREATE TABLE mymodel.networkdep AS

SELECT filename, gid, ST_Value(rast, geom) val

FROM mymodel.deposition, mymodel.network

WHERE ST_Intersects(rast, geom)

ORDER BY gid, rid;


Gathering from the count it looks like all the pixels of those are treated as 
data when some should be treated as no data.  If it's not an error in the 
raster you are using,

it would really help to get output of the polygon and also the other stats I 
mentioned which (st).* outputs  So a query something like:

WITH  foo AS (
  SELECT  mymodel.concentrated.rid,  mymodel.networkpoly.gid,  geom,  
ST_SummaryStats( ST_Clip(rast, geom) ) As st
FROM mymodel.concentrated INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )
   WHERE filename='10_inci.tif'
)
SELECT rid, gid, geom, (st).*
FROM foo;


You might actually want to plot the geom on the map overlaid with the problem 
tiles (or original raster) you classified for that region.

My guess is something went wrong in the ST_Reclass and all the pixels in those 
tiles did not get set to 0 (and thus are treated as data doing a summary stats 
of the problem pre-reclassed tiles and looking at the min,max will give you a 
clue about that).
This second point may be moot if it's just you meant to put deposition instead 
of concentration, but good knowledge to keep in mind for future troubleshooting.

Hope that helps,
Regina



___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Dear Regina,

Wow that makes a tremendous difference.   I used this polygon to extract the 
values using a modified form of the query you suggested in an earlier post and 
the whole thing took under five minutes!

No question postgis is the best way to do this.

I cannot thank you enough – but then again I cannot check the answer it gives 
is right ☺  - but at least the methodology is consistent.

Best wishes

Darrel





From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 28 November 2015 05:04
To: 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

For completeness, your vector query would then look like

CREATE TABLE mymodel.network AS

SELECT rid As gid, ST_Polygon(ST_Reclass(rast,  1, '[0-600):0, 
[600-1]:1','1BB', 0) ) As geom
FROM mymodel.concentrated ;



I think you can put in –number in there like - -1000-600.  There was an issue 
way back with negatives but I thnk that was fixed in PostGIS 2.1 something so 
should work

So what that basically does is create a new raster from original setting of 
pixel values from 0 to < 600 to 0 and from >= 600 to 1 to 1 and making that 
a 1BB (1-bit booleanl raster), and then definiing 0 as no data so that when you 
convert it to a multipolygon you'll only vectorize the 600-1 range.

I think ST_polygon against a 1BB is faster than larger band pixel types.

Hope that helps,
Regina


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: Friday, November 27, 2015 10:04 PM
To: 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,

I think the equivalent in PostGIS terminology would be ST_Reclass - 
http://postgis.net/docs/manual-2.2/RT_ST_Reclass.html

And that's a fairly fast operation as I recall.

Hope that helps,
Regina



From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Friday, November 27, 2015 4:09 PM
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Dear Regina,

Yes, I should have been clearer. To be honest I did not fully understand what I 
was trying to do at the outset ☺

The binary raster to point conversion  was an afterthought as I had already 
been extracting data along short transects using a point shape file.  For the 
record arc has a function in the raster calculator which works as a conditional 
such that
CON(raster, true, false, condition) which in my case was simply CON(FA, 1, 0 
“value > 600”).  This produces a raster which I then had to convert to points 
and eliminate the zeros.  I’m sure there is a way to do this in QGIS but the 
current raster calculator does not allow that directly.  FYI my rasters have 
only one band as they are being used to store numerical model matrix outputs so 
that I can readily visualise them in order to allow me to see structures I 
would not recognise hidden within the 5 million cell datasets.

Anyhow I will experiment with doing more of this in postgis (it would be great 
if I could script an end-to-end solution). Once I am happy with the numerical 
model I will have the model write the data directly into postgis.

I have made some progress this week  thanks to your help. Hopefully I am 
beginning to see how best to use this tool for my intended purpose.

Thanks again

Darrel



From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 27 November 2015 20:02
To: 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,

Oh that's what you are trying to do sorry I didn't recognize that whole CASE 
thing as a binary check operation until you described the purpose.

For the bit operation type stuff it is much faster to define that 0/1 as a 
geometry (which it looks like you've done, but I don't know if you just have 
one pixel cell per or what or details of how you do it in ArcGIS.  I suspect 
that logic can be recreated easily in PostGIS with something like below:

If your network raster is just a set of 0s and 1s (is it a 1BB) or you just 
want to treat 0 as no-date (which is essentially what you case statement was 
trying to do I think) then you could just convert it to a geometry with these 
functions
http://postgis.net/docs/manual-2.2/RT_ST_Polygon.html, 
http://postgis.net/docs/manual-2.2/RT_ST_SetBandNoDataValue.html

and this SQL Statement

CREATE TABLE mymodel.network AS
SELECT rid As gid, ST_Polygon(ST_SetBandNoDataValue(rast,1, 0) ) As geom
FROM mymodel.concentrated ;


Once you have that concentrated as a network channel as a geom

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Oh dear I spoke too soon!

I created the polygon as directed and this is fine when I visualise it in QGIS.

Now when I did the intersects and sum yesterday using the point data for one of 
the rasters I got a total of 0.946094531.

When I do the extract in QGIS using the polygon created in postgis  I get the 
same total.

Alas when I do this (note this is the same raster/polygon intersection):

WITH  foo AS (
  SELECT  mymodel.concentrated.rid,   ST_SummaryStats( ST_Clip(rast, geom) ) As 
st
FROM mymodel.concentrated INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )
   WHERE filename='10_inci.tif'
)
SELECT SUM( (st).sum )
FROM foo;

I get this 21.261754843969

I believe this answer is wrong – so what did I miss?

Thanks

Darrel






From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: 28 November 2015 22:51
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?


Dear Regina,

Wow that makes a tremendous difference.   I used this polygon to extract the 
values using a modified form of the query you suggested in an earlier post and 
the whole thing took under five minutes!

No question postgis is the best way to do this.

I cannot thank you enough – but then again I cannot check the answer it gives 
is right ☺  - but at least the methodology is consistent.

Best wishes

Darrel





From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 28 November 2015 05:04
To: 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

For completeness, your vector query would then look like

CREATE TABLE mymodel.network AS

SELECT rid As gid, ST_Polygon(ST_Reclass(rast,  1, '[0-600):0, 
[600-1]:1','1BB', 0) ) As geom
FROM mymodel.concentrated ;



I think you can put in –number in there like - -1000-600.  There was an issue 
way back with negatives but I thnk that was fixed in PostGIS 2.1 something so 
should work

So what that basically does is create a new raster from original setting of 
pixel values from 0 to < 600 to 0 and from >= 600 to 1 to 1 and making that 
a 1BB (1-bit booleanl raster), and then definiing 0 as no data so that when you 
convert it to a multipolygon you'll only vectorize the 600-1 range.

I think ST_polygon against a 1BB is faster than larger band pixel types.

Hope that helps,
Regina


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: Friday, November 27, 2015 10:04 PM
To: 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,

I think the equivalent in PostGIS terminology would be ST_Reclass - 
http://postgis.net/docs/manual-2.2/RT_ST_Reclass.html

And that's a fairly fast operation as I recall.

Hope that helps,
Regina



From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Friday, November 27, 2015 4:09 PM
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Dear Regina,

Yes, I should have been clearer. To be honest I did not fully understand what I 
was trying to do at the outset ☺

The binary raster to point conversion  was an afterthought as I had already 
been extracting data along short transects using a point shape file.  For the 
record arc has a function in the raster calculator which works as a conditional 
such that
CON(raster, true, false, condition) which in my case was simply CON(FA, 1, 0 
“value > 600”).  This produces a raster which I then had to convert to points 
and eliminate the zeros.  I’m sure there is a way to do this in QGIS but the 
current raster calculator does not allow that directly.  FYI my rasters have 
only one band as they are being used to store numerical model matrix outputs so 
that I can readily visualise them in order to allow me to see structures I 
would not recognise hidden within the 5 million cell datasets.

Anyhow I will experiment with doing more of this in postgis (it would be great 
if I could script an end-to-end solution). Once I am happy with the numerical 
model I will have the model write the data directly into postgis.

I have made some progress this week  thanks to your help. Hopefully I am 
beginning to see how best to use this tool for my intended purpose.

Thanks again

Darrel



From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 27 November 2015 20:02
To: 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,

Oh that's what you a

Re: [postgis-users] Help with SQL query?

[Paragon Corporation]

Darrel,

 

Try not summing.  That might give you more information about  what is wrong.  
So instead do this:

 

WITH  foo AS (

  SELECT  mymodel.concentrated.rid,  mymodel.networkpoly.gid,  ST_SummaryStats( 
ST_Clip(rast, geom) ) As st

FROM mymodel.concentrated INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )

   WHERE filename='10_inci.tif'

)

SELECT rid, gid, (st).*

FROM foo;

 

That might give you multiple records per rid which might be okay since a single 
tile might be clipped by more than one polygon.  At anyrate, you'll get a lot 
more stats and be able to target the tile/poly combo that might be at issue.

 

So your query should have

 

rid, gid, count, sum, mean, stddev, min, max  (for each tile/geom combo – which 
should hopefully give you a clue what went wrong).

 

Hope that helps,

Regina

 

 

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Saturday, November 28, 2015 6:43 PM
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

 

Oh dear I spoke too soon!

 

I created the polygon as directed and this is fine when I visualise it in QGIS. 

 

Now when I did the intersects and sum yesterday using the point data for one of 
the rasters I got a total of 0.946094531.

 

When I do the extract in QGIS using the polygon created in postgis  I get the 
same total.

 

Alas when I do this (note this is the same raster/polygon intersection):

 

WITH  foo AS (

  SELECT  mymodel.concentrated.rid,   ST_SummaryStats( ST_Clip(rast, geom) ) As 
st

FROM mymodel.concentrated INNER JOIN mymodel.networkpoly ON ( 
ST_Intersects(rast,geom) )

   WHERE filename='10_inci.tif'

)

SELECT SUM( (st).sum )

FROM foo;

 

I get this 21.261754843969

 

I believe this answer is wrong – so what did I miss?

 

Thanks

 

Darrel

 

 

 

 

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: 28 November 2015 22:51
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org 
<mailto:postgis-users@lists.osgeo.org> >
Subject: Re: [postgis-users] Help with SQL query?

 

 

Dear Regina,

 

Wow that makes a tremendous difference.   I used this polygon to extract the 
values using a modified form of the query you suggested in an earlier post and 
the whole thing took under five minutes!  

 

No question postgis is the best way to do this.

 

I cannot thank you enough – but then again I cannot check the answer it gives 
is right :)  - but at least the methodology is consistent. 

 

Best wishes

 

Darrel

 

 

 

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 28 November 2015 05:04
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org 
<mailto:postgis-users@lists.osgeo.org> >
Subject: Re: [postgis-users] Help with SQL query?

 

For completeness, your vector query would then look like

 

CREATE TABLE mymodel.network AS

SELECT rid As gid, ST_Polygon(ST_Reclass(rast,  1, '[0-600):0, 
[600-1]:1','1BB', 0) ) As geom

FROM mymodel.concentrated ;

 

 

 

I think you can put in –number in there like - -1000-600.  There was an issue 
way back with negatives but I thnk that was fixed in PostGIS 2.1 something so 
should work

 

So what that basically does is create a new raster from original setting of 
pixel values from 0 to < 600 to 0 and from >= 600 to 1 to 1 and making that 
a 1BB (1-bit booleanl raster), and then definiing 0 as no data so that when you 
convert it to a multipolygon you'll only vectorize the 600-1 range.

 

I think ST_polygon against a 1BB is faster than larger band pixel types.

 

Hope that helps,

Regina

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: Friday, November 27, 2015 10:04 PM
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org 
<mailto:postgis-users@lists.osgeo.org> >
Subject: Re: [postgis-users] Help with SQL query?

 

Darrel,

 

I think the equivalent in PostGIS terminology would be ST_Reclass -  
<http://postgis.net/docs/manual-2.2/RT_ST_Reclass.html> 
http://postgis.net/docs/manual-2.2/RT_ST_Reclass.html

 

And that's a fairly fast operation as I recall.

 

Hope that helps,

Regina

 

 

 

From: postgis-users [ <mailto:postgis-users-boun...@lists.osgeo.org> 
mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of Darrel Maddy
Sent: Friday, November 27, 2015 4:09 PM
To: PostGIS Users Discussion < <mailto:postgis-users@lists.osgeo.org> 
postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

 

Dear Regina,

 

Yes, I should have been clearer. To be honest I did not fully understand what I 
was trying to do at the outset :)

 

The binary raster to point conversion  was an afterthought as

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Dear Regina,

Many thanks for this suggestion.

I ran the query in this form and for one raster it takes 3372s (~55 mins).  I 
guess that is what I had anticipated from the single tile  exercise I ran with 
the alternate algorithm.

This still seems a little too long however and so I have started to explore 
ways to improve upon this by pre-processing some of the data.  The cells of 
interest represent the ‘main channels’ in a drainage network.  Although in 
future this network may change position from one output timestep to another 
(they are actually 1000l iterations of the model apart) in this particular 
variant the position of these cells is static.  With that in mind I decided to 
create a binary raster where the 1’s represent the channel cells (I had to do 
this in arcgis as QGIS does not appear to have a Con function in the raster 
calculator!)  and then exported this as a point layer. I then deleted the 
points coded zero and saved the shp file in QGIS.  I imported the ntwork shp 
file into postgis (there are about 15500 15 points) and I am now running:

CREATE TABLE mymodel.networkdep AS
SELECT filename, gid, ST_Value(rast, geom) val
FROM mymodel.deposition, mymodel.network
WHERE ST_Intersects(rast, geom)
ORDER BY gid, rid;

I will sum by raster (i.e. filename) using the new table but will settle for 
just having the relevant data for now. This took  2057s(~35 minutes!) to 
complete!

If this is the best way to do it I will explore the OGR library and try and 
hardcode the network point output directly into the model or, more 
realistically, write a short routine to extract this automatically from the 
flow accumulation output rasters without recourse to a GIS.

I am learning a lot through this exercise, so thanks once again to all of you 
who have made suggestions, they are very much appreciated.

Best wishes

Darrel








From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 26 November 2015 05:08
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

That timing seems much slower than I recall.

FWIW expression based mapalgebra as I recall is slower than using the call back 
function approach.  So you could try wrapping your CASE in a call back function.

However I think something else might be going on here and postgres might be 
repeating work.  I forgot under what conditions it decides to reevaluate a 
function call, I just remember being really surprised by it.

To avoid that, you can try using a CTE, also you don't need that ST_Union call 
which for larger number of rasters is expensive, and you might even generate a 
raster that is too big to compute.

I'm also guessing your rasts are all tiled the same, so you really don't need 
ST_Intersects, just use the same box operator

So try this:

WITH  foo AS (
  SELECT ST_SummaryStats( ST_MapAlgebra(deposition.rast, concentrated.rast, 
'CASE WHEN [rast2] > 0.0 THEN [rast1] ELSE NULL END' ) ) As st
FROM mymodel.deposition INNER JOIN mymodel.concentrated ON ( 
deposition.rast  ~=  concentrated.rast )
WHERE deposition.rid=1

)
SELECT SUM( (st).sum )
FROM foo;


Hope that helps,
Regina
http://www.postgis.us
http://postgis.net


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Wednesday, November 25, 2015 5:06 PM
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>; Brent 
Wood <pcr...@yahoo.com<mailto:pcr...@yahoo.com>>
Subject: Re: [postgis-users] Help with SQL query?

Dear Brent,

I must confess that my attempts to do this are so far proving very unsuccessful

If  I run the following query:

SELECT  (ST_SummaryStats(ST_Union(rast))).sum AS sum
FROM  (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'CASE WHEN 
[rast2] > 0.0 THEN [rast1] ELSE NULL END' ) As rast
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_Intersects(deposition.rast, concentrated.rast) AND 
deposition.rid=1 ) foo ;

It takes around 30 seconds to complete as I assume it is only looking at one 
tile(they are 256x256 pixels) i.e. rid 1. It is not easy to check the sum – for 
that I need one complete raster.

For the record this was marginally faster than
SELECT (ST_SummaryStats(ST_Union(rast))).sum AS sum
FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'CASE WHEN 
[rast2] > 0.0 THEN [rast1] ELSE NULL END' ) As rast
FROM mymodel.deposition, mymodel.concentrated
WHERE mymodel.deposition.filename='10_depo.tif' AND 
ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast) ) foo ;
Even after I built indexes for the clauses after the WHERE.

Now there are 144 tiles in each of the rasters I want to p

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Dear Regina,

Yes, I should have been clearer. To be honest I did not fully understand what I 
was trying to do at the outset ☺

The binary raster to point conversion  was an afterthought as I had already 
been extracting data along short transects using a point shape file.  For the 
record arc has a function in the raster calculator which works as a conditional 
such that
CON(raster, true, false, condition) which in my case was simply CON(FA, 1, 0 
“value > 600”).  This produces a raster which I then had to convert to points 
and eliminate the zeros.  I’m sure there is a way to do this in QGIS but the 
current raster calculator does not allow that directly.  FYI my rasters have 
only one band as they are being used to store numerical model matrix outputs so 
that I can readily visualise them in order to allow me to see structures I 
would not recognise hidden within the 5 million cell datasets.

Anyhow I will experiment with doing more of this in postgis (it would be great 
if I could script an end-to-end solution). Once I am happy with the numerical 
model I will have the model write the data directly into postgis.

I have made some progress this week  thanks to your help. Hopefully I am 
beginning to see how best to use this tool for my intended purpose.

Thanks again

Darrel



From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 27 November 2015 20:02
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,

Oh that's what you are trying to do sorry I didn't recognize that whole CASE 
thing as a binary check operation until you described the purpose.

For the bit operation type stuff it is much faster to define that 0/1 as a 
geometry (which it looks like you've done, but I don't know if you just have 
one pixel cell per or what or details of how you do it in ArcGIS.  I suspect 
that logic can be recreated easily in PostGIS with something like below:

If your network raster is just a set of 0s and 1s (is it a 1BB) or you just 
want to treat 0 as no-date (which is essentially what you case statement was 
trying to do I think) then you could just convert it to a geometry with these 
functions
http://postgis.net/docs/manual-2.2/RT_ST_Polygon.html, 
http://postgis.net/docs/manual-2.2/RT_ST_SetBandNoDataValue.html

and this SQL Statement

CREATE TABLE mymodel.network AS
SELECT rid As gid, ST_Polygon(ST_SetBandNoDataValue(rast,1, 0) ) As geom
FROM mymodel.concentrated ;


Once you have that concentrated as a network channel as a geometry, then you 
can use ST_Clip and that should be pretty fast and give you the same results.
http://postgis.net/docs/manual-2.2/RT_ST_Clip.html


So your query would look something like

WITH  foo AS (
  SELECT  mymodel.deposition.rid,   ST_SummaryStats( ST_Clip(rast, geom) ) As st
FROM mymodel.deposition INNER JOIN mymodel.network ON ( 
ST_Intersects(rast,geom) )

)
SELECT SUM( (st).sum )
FROM foo;

Hope that helps,
Regina
http://www.postgis.us
http://postgis.net


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Friday, November 27, 2015 12:42 PM
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Dear Regina,

Many thanks for this suggestion.

I ran the query in this form and for one raster it takes 3372s (~55 mins).  I 
guess that is what I had anticipated from the single tile  exercise I ran with 
the alternate algorithm.

This still seems a little too long however and so I have started to explore 
ways to improve upon this by pre-processing some of the data.  The cells of 
interest represent the ‘main channels’ in a drainage network.  Although in 
future this network may change position from one output timestep to another 
(they are actually 1000l iterations of the model apart) in this particular 
variant the position of these cells is static.  With that in mind I decided to 
create a binary raster where the 1’s represent the channel cells (I had to do 
this in arcgis as QGIS does not appear to have a Con function in the raster 
calculator!)  and then exported this as a point layer. I then deleted the 
points coded zero and saved the shp file in QGIS.  I imported the ntwork shp 
file into postgis (there are about 15500 15 points) and I am now running:

CREATE TABLE mymodel.networkdep AS
SELECT filename, gid, ST_Value(rast, geom) val
FROM mymodel.deposition, mymodel.network
WHERE ST_Intersects(rast, geom)
ORDER BY gid, rid;

I will sum by raster (i.e. filename) using the new table but will settle for 
just having the relevant data for now. This took  2057s(~35 minutes!) to 
complete!

If this is the best way to do it I will explore the OGR library and try and 
hardcode the network point output directly into the model or, more 
realistically, write a short routine to extra

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Dear Brent,

Yes that works for me too !

Sometimes I look for what is familiar and overlook the obvious – a little 
embarrassing, but I know now .

Many thanks

Darrel

From: Brent Wood [mailto:pcr...@yahoo.com]
Sent: 27 November 2015 21:53
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org>; Darrel Maddy 
<darrel.ma...@newcastle.ac.uk>
Subject: Re: [postgis-users] Help with SQL query?


Darrel

You can use the logical operators in the QGIS raster calculator as a binary 
conditional :

eg: raster a is used to generate a new raster with default 0, but pixel value=1 
 where pixels in a have a value >300...   use the following operation

(a > 300.0)

Save the result as a grid rather than a image. It works for me
Brent

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Paragon Corporation]

Darrel,

 

Oh that's what you are trying to do sorry I didn't recognize that whole CASE 
thing as a binary check operation until you described the purpose.

 

For the bit operation type stuff it is much faster to define that 0/1 as a 
geometry (which it looks like you've done, but I don't know if you just have 
one pixel cell per or what or details of how you do it in ArcGIS.  I suspect 
that logic can be recreated easily in PostGIS with something like below:

 

If your network raster is just a set of 0s and 1s (is it a 1BB) or you just 
want to treat 0 as no-date (which is essentially what you case statement was 
trying to do I think) then you could just convert it to a geometry with these 
functions

http://postgis.net/docs/manual-2.2/RT_ST_Polygon.html, 
http://postgis.net/docs/manual-2.2/RT_ST_SetBandNoDataValue.html 

 

and this SQL Statement

 

CREATE TABLE mymodel.network AS

SELECT rid As gid, ST_Polygon(ST_SetBandNoDataValue(rast,1, 0) ) As geom

FROM mymodel.concentrated ;

 

 

Once you have that concentrated as a network channel as a geometry, then you 
can use ST_Clip and that should be pretty fast and give you the same results.

http://postgis.net/docs/manual-2.2/RT_ST_Clip.html

 

 

So your query would look something like

 

WITH  foo AS (

  SELECT  mymodel.deposition.rid,   ST_SummaryStats( ST_Clip(rast, geom) ) As st

FROM mymodel.deposition INNER JOIN mymodel.network ON ( 
ST_Intersects(rast,geom) )

   

)

SELECT SUM( (st).sum )

FROM foo;

 

Hope that helps,

Regina

http://www.postgis.us

http://postgis.net

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Friday, November 27, 2015 12:42 PM
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

 

Dear Regina,

 

Many thanks for this suggestion.

 

I ran the query in this form and for one raster it takes 3372s (~55 mins).  I 
guess that is what I had anticipated from the single tile  exercise I ran with 
the alternate algorithm.

 

This still seems a little too long however and so I have started to explore 
ways to improve upon this by pre-processing some of the data.  The cells of 
interest represent the ‘main channels’ in a drainage network.  Although in 
future this network may change position from one output timestep to another 
(they are actually 1000l iterations of the model apart) in this particular 
variant the position of these cells is static.  With that in mind I decided to 
create a binary raster where the 1’s represent the channel cells (I had to do 
this in arcgis as QGIS does not appear to have a Con function in the raster 
calculator!)  and then exported this as a point layer. I then deleted the 
points coded zero and saved the shp file in QGIS.  I imported the ntwork shp 
file into postgis (there are about 15500 15 points) and I am now running:

 

CREATE TABLE mymodel.networkdep AS

SELECT filename, gid, ST_Value(rast, geom) val

FROM mymodel.deposition, mymodel.network

WHERE ST_Intersects(rast, geom) 

ORDER BY gid, rid;

 

I will sum by raster (i.e. filename) using the new table but will settle for 
just having the relevant data for now. This took  2057s(~35 minutes!) to 
complete! 

 

If this is the best way to do it I will explore the OGR library and try and 
hardcode the network point output directly into the model or, more 
realistically, write a short routine to extract this automatically from the 
flow accumulation output rasters without recourse to a GIS.

 

I am learning a lot through this exercise, so thanks once again to all of you 
who have made suggestions, they are very much appreciated.

 

Best wishes

 

Darrel 

 

 

 

 

 

 

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 26 November 2015 05:08
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org 
<mailto:postgis-users@lists.osgeo.org> >
Subject: Re: [postgis-users] Help with SQL query?

 

That timing seems much slower than I recall.

 

FWIW expression based mapalgebra as I recall is slower than using the call back 
function approach.  So you could try wrapping your CASE in a call back function.

 

However I think something else might be going on here and postgres might be 
repeating work.  I forgot under what conditions it decides to reevaluate a 
function call, I just remember being really surprised by it.

 

To avoid that, you can try using a CTE, also you don't need that ST_Union call 
which for larger number of rasters is expensive, and you might even generate a 
raster that is too big to compute.

 

I'm also guessing your rasts are all tiled the same, so you really don't need 
ST_Intersects, just use the same box operator

 

So try this:

 

WITH  foo AS (

  SELECT ST_SummaryStats( ST_MapAlgebra(deposition.rast, concentrated.rast, 
'CASE WHEN [rast2] > 0.0 THEN [rast1] ELS

Re: [postgis-users] Help with SQL query?

[Paragon Corporation]

Darrel,

 

I think the equivalent in PostGIS terminology would be ST_Reclass - 
http://postgis.net/docs/manual-2.2/RT_ST_Reclass.html

 

And that's a fairly fast operation as I recall.

 

Hope that helps,

Regina

 

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Friday, November 27, 2015 4:09 PM
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

 

Dear Regina,

 

Yes, I should have been clearer. To be honest I did not fully understand what I 
was trying to do at the outset :)

 

The binary raster to point conversion  was an afterthought as I had already 
been extracting data along short transects using a point shape file.  For the 
record arc has a function in the raster calculator which works as a conditional 
such that 
CON(raster, true, false, condition) which in my case was simply CON(FA, 1, 0 
“value > 600”).  This produces a raster which I then had to convert to points 
and eliminate the zeros.  I’m sure there is a way to do this in QGIS but the 
current raster calculator does not allow that directly.  FYI my rasters have 
only one band as they are being used to store numerical model matrix outputs so 
that I can readily visualise them in order to allow me to see structures I 
would not recognise hidden within the 5 million cell datasets.

 

Anyhow I will experiment with doing more of this in postgis (it would be great 
if I could script an end-to-end solution). Once I am happy with the numerical 
model I will have the model write the data directly into postgis. 

 

I have made some progress this week  thanks to your help. Hopefully I am 
beginning to see how best to use this tool for my intended purpose.

 

Thanks again

 

Darrel

 

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: 27 November 2015 20:02
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org 
<mailto:postgis-users@lists.osgeo.org> >
Subject: Re: [postgis-users] Help with SQL query?

 

Darrel,

 

Oh that's what you are trying to do sorry I didn't recognize that whole CASE 
thing as a binary check operation until you described the purpose.

 

For the bit operation type stuff it is much faster to define that 0/1 as a 
geometry (which it looks like you've done, but I don't know if you just have 
one pixel cell per or what or details of how you do it in ArcGIS.  I suspect 
that logic can be recreated easily in PostGIS with something like below:

 

If your network raster is just a set of 0s and 1s (is it a 1BB) or you just 
want to treat 0 as no-date (which is essentially what you case statement was 
trying to do I think) then you could just convert it to a geometry with these 
functions

http://postgis.net/docs/manual-2.2/RT_ST_Polygon.html, 
http://postgis.net/docs/manual-2.2/RT_ST_SetBandNoDataValue.html 

 

and this SQL Statement

 

CREATE TABLE mymodel.network AS

SELECT rid As gid, ST_Polygon(ST_SetBandNoDataValue(rast,1, 0) ) As geom

FROM mymodel.concentrated ;

 

 

Once you have that concentrated as a network channel as a geometry, then you 
can use ST_Clip and that should be pretty fast and give you the same results.

http://postgis.net/docs/manual-2.2/RT_ST_Clip.html

 

 

So your query would look something like

 

WITH  foo AS (

  SELECT  mymodel.deposition.rid,   ST_SummaryStats( ST_Clip(rast, geom) ) As st

FROM mymodel.deposition INNER JOIN mymodel.network ON ( 
ST_Intersects(rast,geom) )

   

)

SELECT SUM( (st).sum )

FROM foo;

 

Hope that helps,

Regina

http://www.postgis.us

http://postgis.net

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Friday, November 27, 2015 12:42 PM
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org 
<mailto:postgis-users@lists.osgeo.org> >
Subject: Re: [postgis-users] Help with SQL query?

 

Dear Regina,

 

Many thanks for this suggestion.

 

I ran the query in this form and for one raster it takes 3372s (~55 mins).  I 
guess that is what I had anticipated from the single tile  exercise I ran with 
the alternate algorithm.

 

This still seems a little too long however and so I have started to explore 
ways to improve upon this by pre-processing some of the data.  The cells of 
interest represent the ‘main channels’ in a drainage network.  Although in 
future this network may change position from one output timestep to another 
(they are actually 1000l iterations of the model apart) in this particular 
variant the position of these cells is static.  With that in mind I decided to 
create a binary raster where the 1’s represent the channel cells (I had to do 
this in arcgis as QGIS does not appear to have a Con function in the raster 
calculator!)  and then exported this as a point layer. I then deleted the 
points coded zero and saved the shp file in QGIS.

Re: [postgis-users] Help with SQL query?

[Paragon Corporation]

For completeness, your vector query would then look like

 

CREATE TABLE mymodel.network AS

SELECT rid As gid, ST_Polygon(ST_Reclass(rast,  1, '[0-600):0, 
[600-1]:1','1BB', 0) ) As geom

FROM mymodel.concentrated ;

 

 

 

I think you can put in –number in there like - -1000-600.  There was an issue 
way back with negatives but I thnk that was fixed in PostGIS 2.1 something so 
should work

 

So what that basically does is create a new raster from original setting of 
pixel values from 0 to < 600 to 0 and from >= 600 to 1 to 1 and making that 
a 1BB (1-bit booleanl raster), and then definiing 0 as no data so that when you 
convert it to a multipolygon you'll only vectorize the 600-1 range.

 

I think ST_polygon against a 1BB is faster than larger band pixel types.

 

Hope that helps,

Regina

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Paragon Corporation
Sent: Friday, November 27, 2015 10:04 PM
To: 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

 

Darrel,

 

I think the equivalent in PostGIS terminology would be ST_Reclass -  
<http://postgis.net/docs/manual-2.2/RT_ST_Reclass.html> 
http://postgis.net/docs/manual-2.2/RT_ST_Reclass.html

 

And that's a fairly fast operation as I recall.

 

Hope that helps,

Regina

 

 

 

From: postgis-users [ <mailto:postgis-users-boun...@lists.osgeo.org> 
mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of Darrel Maddy
Sent: Friday, November 27, 2015 4:09 PM
To: PostGIS Users Discussion < <mailto:postgis-users@lists.osgeo.org> 
postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

 

Dear Regina,

 

Yes, I should have been clearer. To be honest I did not fully understand what I 
was trying to do at the outset :)

 

The binary raster to point conversion  was an afterthought as I had already 
been extracting data along short transects using a point shape file.  For the 
record arc has a function in the raster calculator which works as a conditional 
such that 
CON(raster, true, false, condition) which in my case was simply CON(FA, 1, 0 
“value > 600”).  This produces a raster which I then had to convert to points 
and eliminate the zeros.  I’m sure there is a way to do this in QGIS but the 
current raster calculator does not allow that directly.  FYI my rasters have 
only one band as they are being used to store numerical model matrix outputs so 
that I can readily visualise them in order to allow me to see structures I 
would not recognise hidden within the 5 million cell datasets.

 

Anyhow I will experiment with doing more of this in postgis (it would be great 
if I could script an end-to-end solution). Once I am happy with the numerical 
model I will have the model write the data directly into postgis. 

 

I have made some progress this week  thanks to your help. Hopefully I am 
beginning to see how best to use this tool for my intended purpose.

 

Thanks again

 

Darrel

 

 

 

From: postgis-users [ <mailto:postgis-users-boun...@lists.osgeo.org> 
mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of Paragon Corporation
Sent: 27 November 2015 20:02
To: 'PostGIS Users Discussion' < <mailto:postgis-users@lists.osgeo.org> 
postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

 

Darrel,

 

Oh that's what you are trying to do sorry I didn't recognize that whole CASE 
thing as a binary check operation until you described the purpose.

 

For the bit operation type stuff it is much faster to define that 0/1 as a 
geometry (which it looks like you've done, but I don't know if you just have 
one pixel cell per or what or details of how you do it in ArcGIS.  I suspect 
that logic can be recreated easily in PostGIS with something like below:

 

If your network raster is just a set of 0s and 1s (is it a 1BB) or you just 
want to treat 0 as no-date (which is essentially what you case statement was 
trying to do I think) then you could just convert it to a geometry with these 
functions

 <http://postgis.net/docs/manual-2.2/RT_ST_Polygon.html> 
http://postgis.net/docs/manual-2.2/RT_ST_Polygon.html,  
<http://postgis.net/docs/manual-2.2/RT_ST_SetBandNoDataValue.html> 
http://postgis.net/docs/manual-2.2/RT_ST_SetBandNoDataValue.html 

 

and this SQL Statement

 

CREATE TABLE mymodel.network AS

SELECT rid As gid, ST_Polygon(ST_SetBandNoDataValue(rast,1, 0) ) As geom

FROM mymodel.concentrated ;

 

 

Once you have that concentrated as a network channel as a geometry, then you 
can use ST_Clip and that should be pretty fast and give you the same results.

 <http://postgis.net/docs/manual-2.2/RT_ST_Clip.html> 
http://postgis.net/docs/manual-2.2/RT_ST_Clip.html

 

 

So your query would look something like

 

WITH  foo AS (

  SELECT  mymodel.deposition.rid,   ST_SummaryStats( ST_Clip(rast, ge

Re: [postgis-users] Help with SQL query?

[Paragon Corporation]

That timing seems much slower than I recall.

 

FWIW expression based mapalgebra as I recall is slower than using the call back 
function approach.  So you could try wrapping your CASE in a call back function.

 

However I think something else might be going on here and postgres might be 
repeating work.  I forgot under what conditions it decides to reevaluate a 
function call, I just remember being really surprised by it.

 

To avoid that, you can try using a CTE, also you don't need that ST_Union call 
which for larger number of rasters is expensive, and you might even generate a 
raster that is too big to compute.

 

I'm also guessing your rasts are all tiled the same, so you really don't need 
ST_Intersects, just use the same box operator

 

So try this:

 

WITH  foo AS (

  SELECT ST_SummaryStats( ST_MapAlgebra(deposition.rast, concentrated.rast, 
'CASE WHEN [rast2] > 0.0 THEN [rast1] ELSE NULL END' ) ) As st

FROM mymodel.deposition INNER JOIN mymodel.concentrated ON ( 
deposition.rast  ~=  concentrated.rast )

WHERE deposition.rid=1

 

)

SELECT SUM( (st).sum )

FROM foo;

 

 

Hope that helps,

Regina

http://www.postgis.us

http://postgis.net

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Wednesday, November 25, 2015 5:06 PM
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org>; Brent Wood 
<pcr...@yahoo.com>
Subject: Re: [postgis-users] Help with SQL query?

 

Dear Brent,

 

I must confess that my attempts to do this are so far proving very unsuccessful

 

If  I run the following query:

 

SELECT  (ST_SummaryStats(ST_Union(rast))).sum AS sum

FROM  (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'CASE WHEN 
[rast2] > 0.0 THEN [rast1] ELSE NULL END' ) As rast

FROM mymodel.deposition, mymodel.concentrated

WHERE ST_Intersects(deposition.rast, concentrated.rast) AND 
deposition.rid=1 ) foo ;

 

It takes around 30 seconds to complete as I assume it is only looking at one 
tile(they are 256x256 pixels) i.e. rid 1. It is not easy to check the sum – for 
that I need one complete raster.

 

For the record this was marginally faster than

SELECT (ST_SummaryStats(ST_Union(rast))).sum AS sum

FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'CASE WHEN 
[rast2] > 0.0 THEN [rast1] ELSE NULL END' ) As rast

FROM mymodel.deposition, mymodel.concentrated

WHERE mymodel.deposition.filename='10_depo.tif' AND 
ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND 

 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast) ) foo ;

Even after I built indexes for the clauses after the WHERE.

 

Now there are 144 tiles in each of the rasters I want to perform this operation 
on.  Logic would therefore suggest this should take ~4500s

 

However when I perform the following query

 

SELECT  (ST_SummaryStats(ST_Union(rast))).sum AS sum

FROM  (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'CASE WHEN 
[rast2] > 0.0 THEN [rast1] ELSE NULL END' ) As rast

FROM mymodel.deposition, mymodel.concentrated

WHERE ST_Intersects(deposition.rast, concentrated.rast) AND 
deposition.filename='10_depo.tif' ) foo ;

 

The query is still running after 18000s!  I must therefore assume I have done 
something wrong but as you may have guessed the answer eludes me.

 

Any further suggestions would be welcome but I will continue to try and find a 
solution as I have 135 rasters to perform this operations on now and 
potentially many thousands more in the future.

 

Darrel

 

.

 

I

 

 

 

 

 

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: 24 November 2015 19:52
To: Brent Wood <pcr...@yahoo.com <mailto:pcr...@yahoo.com> >; 
postgis-users@lists.osgeo.org <mailto:postgis-users@lists.osgeo.org> 
Subject: Re: [postgis-users] Help with SQL query?

 

Dear Brent,

 

Many thanks. The data are tiled (256x256) hence the large number of rows from 
the original 135 tifs. I did not build any indexes however, so I will do some 
reading and see how best to approach that (the threads you listed look useful 
so thanks for that).

 

I will run some additional mini queries limited to just one comparison and 
check using QGIS as you suggest – I probably should have done that first!

 

My workstation has 64GB Ram and I would be surprised if it was significantly 
caching to disk. I also have a hexacore intel extreme processor so I would not 
expect this to be hardware limited. I must confess I expected it to finish 
within a couple of hours.

 

Anyhow very many thanks. I will continue to explore and report back hopefully 
with positive news.

 

Darrel

 

 

From: Brent Wood [mailto:pcr...@yahoo.com] 
Sent: 24 November, 2015 7:36 PM
To: Darrel Maddy <da

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Dear Brent,

Many thanks. The data are tiled (256x256) hence the large number of rows from 
the original 135 tifs. I did not build any indexes however, so I will do some 
reading and see how best to approach that (the threads you listed look useful 
so thanks for that).

I will run some additional mini queries limited to just one comparison and 
check using QGIS as you suggest – I probably should have done that first!

My workstation has 64GB Ram and I would be surprised if it was significantly 
caching to disk. I also have a hexacore intel extreme processor so I would not 
expect this to be hardware limited. I must confess I expected it to finish 
within a couple of hours.

Anyhow very many thanks. I will continue to explore and report back hopefully 
with positive news.

Darrel


From: Brent Wood [mailto:pcr...@yahoo.com]
Sent: 24 November, 2015 7:36 PM
To: Darrel Maddy <darrel.ma...@newcastle.ac.uk>; postgis-users@lists.osgeo.org
Subject: Re: [postgis-users] Help with SQL query?

Indexing can improve performance by 100s of x, without them things can be slow. 
Also, did you tile the images when you imported them? If not, then each 
iteration is working through all the pixels in the image, rather than a small 
subset. Essentially with tiles, you have a deep (long) table rather than a wide 
one. RDBMSs work better with lots of small records than a few wide ones, 
especially when indexes are used.

This might help:
http://gis.stackexchange.com/questions/43053/how-to-speed-up-queries-for-raster-databases

and see the raster tutorial they mention for the SRTM data, as to how that is 
loaded into Postgis:
https://trac.osgeo.org/postgis/wiki/WKTRasterTutorial01

To test the logic (the syntax is correct or it wouldn't be working) you could 
add to the "where" clause an extra filter so that only a small subset of the 
entire dataset is included (like just one QGIS operation) then compare this 
with the QGIS result.

That would be much faster that testing on the entire dataset. Once you know it 
is correct for the test case(s), then you can run it on the complete set.

Note that some queries can build up large in-memory objects, so make sure your 
system is not swapping to disk, as that will also slow things down (hugely).

Cheers

Brent
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Dear Roxanne,

Many thanks. I did pickup that issue shortly after my last post and it cured 
the first problem.

There is still an SQL issue with the expression however.

I tried ‘IF [rast2]>0.0 THEN [rast1] ELSE NULL ENDIF’

This produces an error at $2. I can remove that problem by doing this
‘IF ([rast2]>0.0) THEN [rast1] ELSE NULL ENDIF’
But then I get the error at THEN

Curiously ‘[rast2]+[rast1]’  works as intended so somehow the conditional is 
being specified incorrectly.

Unfortunately there is nothing in the documentation which helps me with this 
variant.

I will keep trying.

Darrel

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Roxanne Reid-Bennett
Sent: 24 November, 2015 12:46 AM
To: postgis-users@lists.osgeo.org
Subject: Re: [postgis-users] Help with SQL query?

On 11/23/2015 12:41 PM, Darrel Maddy wrote:
Dear Pierre,

I was not looking for a total solution and I am grateful for the suggestion.

Although it may not look like it, I did consult the documentation and also I 
have Regina’s book beside me.  Unfortunately, for me at least, both documents 
assume some knowledge of SQL – which I do not have.  I am also trying to do 
this simultaneously with a large number of other things that are new to me.  I 
do not find the errors reported at all informative and consider the query I am 
trying to perform to be relatively trivial and hence I had hoped the structure 
of the query might have been more intuitive.  For others it may be.

FWIW - I don't play with rasters, but this appears to be a pure SQL thing... 
add "as rast" like below and try again.


SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF concentrated > 6 
THEN deposition ELSE NULL ENDIF ' )
as rast
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast)

Roxanne



I will try not to bother you again.

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 20:30
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org><mailto:postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

The expression has to stay as it was: 'IF [rast2] > threshold THEN [rast1] ELSE 
NULL ENDIF '

Just replace the threshold value as you did.

Do not try to replace the [rast2] and [rast1]. They refer to the first and 
second raster pixel values. Read the ST_Mapalgebra doc…

Don’t expect our suggestions to work blindly. I did not test this query. I’m 
not in your context. I expect you read the doc about all the mentioned 
functions and adjust for your specific context. I said “your query should “look 
like” this”…

Pierre


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 2:47 PM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] Help with SQL query?

OK I spoke too soon.

I tried this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF concentrated 
> 6 THEN deposition ELSE NULL ENDIF ' )
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast) ) foo;

And all I get is rast does not exist.

I’m afraid the penny has not dropped yet ☹

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: 23 November 2015 16:30
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Dear Pierre and Rasmus,

Many thanks for trying to help.

Rasmus: I am aware of GMT but I was looking for a solution in postgis so that I 
can keep all of the data extraction in one place.

Pierre: That is exactly what I was looking for and very many thanks for 
including the explanation.  I am a little overwhelmed with the number of 
functions offered in postgis. It is certainly a remarkable tool.  Watching the 
queries plough through my datasets is a pleasure – albeit the results do not 
always please me ☺

Hopefully I can put this to work later tonight.

Best wishes

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 16:17
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,


1)  Create a new raster selecting the right pixels with 
ST_MapAlgebra(raster, raster)

2)  Make sure only 

Re: [postgis-users] Help with SQL query?

[Birgit Laggner]

Hi Darrel,

my PostGIS version is too old for testing, but if I read the 
documentation right, then your expression has to be SQL. And IF THEN 
ELSE etc. is not SQL as far as I know - SQL has CASE WHEN.


So, I would assume, you would need to write your expression like this:

'CASE WHEN [rast2] > 0.0 THEN [rast1] ELSE NULL END'

I am curious if this helps with your error...

Regards,

Birgit


Am 24.11.2015 um 10:09 schrieb Darrel Maddy:


Dear Roxanne,

Many thanks. I did pickup that issue shortly after my last post and it 
cured the first problem.


There is still an SQL issue with the expression however.

I tried ‘IF [rast2]>0.0 THEN [rast1] ELSE NULL ENDIF’

This produces an error at $2. I can remove that problem by doing this

‘IF ([rast2]>0.0) THEN [rast1] ELSE NULL ENDIF’

But then I get the error at THEN

Curiously ‘[rast2]+[rast1]’  works as intended so somehow the 
conditional is being specified incorrectly.


Unfortunately there is nothing in the documentation which helps me 
with this variant.


I will keep trying.

Darrel

*From:*postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] 
*On Behalf Of *Roxanne Reid-Bennett

*Sent:* 24 November, 2015 12:46 AM
*To:* postgis-users@lists.osgeo.org
*Subject:* Re: [postgis-users] Help with SQL query?

On 11/23/2015 12:41 PM, Darrel Maddy wrote:

Dear Pierre,

I was not looking for a total solution and I am grateful for the
suggestion.

Although it may not look like it, I did consult the documentation
and also I have Regina’s book beside me.  Unfortunately, for me at
least, both documents assume some knowledge of SQL – which I do
not have.  I am also trying to do this simultaneously with a large
number of other things that are new to me.  I do not find the
errors reported at all informative and consider the query I am
trying to perform to be relatively trivial and hence I had hoped
the structure of the query might have been more intuitive.  For
others it may be.


FWIW - I don't play with rasters, but this appears to be a pure SQL 
thing... add "as rast" like below and try again.



SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF 
concentrated > 6 THEN deposition ELSE NULL ENDIF ' )

as rast

FROM mymodel.deposition, mymodel.concentrated

WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND


 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast)


Roxanne



I will try not to bother you again.

Darrel

*From:*postgis-users
[mailto:postgis-users-boun...@lists.osgeo.org] *On Behalf Of
*Pierre Racine
*Sent:* 23 November 2015 20:30
*To:* PostGIS Users Discussion <postgis-users@lists.osgeo.org>
<mailto:postgis-users@lists.osgeo.org>
    *Subject:* Re: [postgis-users] Help with SQL query?

The expression has to stay as it was: 'IF [rast2] > threshold THEN
[rast1] ELSE NULL ENDIF '

Just replace the threshold value as you did.

Do not try to replace the [rast2] and [rast1]. They refer to the
first and second raster pixel values. Read the ST_Mapalgebra doc…

Don’t expect our suggestions to work blindly. I did not test this
query. I’m not in your context. I expect you read the doc about
all the mentioned functions and adjust for your specific context.
I said “your query should “look like” this”…

Pierre

*From:*postgis-users
[mailto:postgis-users-boun...@lists.osgeo.org] *On Behalf Of
*Darrel Maddy
*Sent:* Monday, November 23, 2015 2:47 PM
*To:* PostGIS Users Discussion
    *Subject:* Re: [postgis-users] Help with SQL query?

OK I spoke too soon.

I tried this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum

FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF
concentrated > 6 THEN deposition ELSE NULL ENDIF ' )

FROM mymodel.deposition, mymodel.concentrated

WHERE ST_UpperleftX(mymodel.deposition.rast) =
ST_UpperleftX(mymodel.concentrated.rast) AND

 ST_UpperleftY(mymodel.deposition.rast) =
ST_UpperleftY(mymodel.deposition.rast) ) foo;

And all I get is rast does not exist.

I’m afraid the penny has not dropped yet L

Darrel

*From:*postgis-users
[mailto:postgis-users-boun...@lists.osgeo.org] *On Behalf Of
*Darrel Maddy
*Sent:* 23 November 2015 16:30
*To:* PostGIS Users Discussion <postgis-users@lists.osgeo.org
<mailto:postgis-users@lists.osgeo.org>>
    *Subject:* Re: [postgis-users] Help with SQL query?

Dear Pierre and Rasmus,

Many thanks for trying to help.

Rasmus: I am aware of GMT but I was looking for a solution in
postgis so that I can keep all of the data extraction in one place.

Pierre: That is exactly what I was looking for and very many
thanks for including the explanation.  I am a little overwhelmed

Re: [postgis-users] Help with SQL query?

[Birgit Laggner]

Thanks for your feedback, Darrel! I hope, your query results will look fine.

Regards,

Birgit


Am 24.11.2015 um 14:00 schrieb Darrel Maddy:


Dear Birgit,

Very many thanks. The query is now running – so it is no longer giving 
an error. I will not know for a few hours whether it is calculating 
the necessary values but it seems progress is being made.


I have a book coming on SQL – I clearly need to do my homework more 
thoroughly J


I will confirm when it stops!

Best wishes

Darrel

*From:*postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] 
*On Behalf Of *Birgit Laggner

*Sent:* 24 November 2015 12:19
*To:* postgis-users@lists.osgeo.org
*Subject:* Re: [postgis-users] Help with SQL query?

Hi Darrel,

my PostGIS version is too old for testing, but if I read the 
documentation right, then your expression has to be SQL. And IF THEN 
ELSE etc. is not SQL as far as I know - SQL has CASE WHEN.


So, I would assume, you would need to write your expression like this:

'CASE WHEN [rast2] > 0.0 THEN [rast1] ELSE NULL END'

I am curious if this helps with your error...

Regards,

Birgit

Am 24.11.2015 um 10:09 schrieb Darrel Maddy:

Dear Roxanne,

Many thanks. I did pickup that issue shortly after my last post
and it cured the first problem.

There is still an SQL issue with the expression however.

I tried ‘IF [rast2]>0.0 THEN [rast1] ELSE NULL ENDIF’

This produces an error at $2. I can remove that problem by doing this

‘IF ([rast2]>0.0) THEN [rast1] ELSE NULL ENDIF’

But then I get the error at THEN

Curiously ‘[rast2]+[rast1]’  works as intended so somehow the
conditional is being specified incorrectly.

Unfortunately there is nothing in the documentation which helps me
with this variant.

I will keep trying.

Darrel

*From:*postgis-users
[mailto:postgis-users-boun...@lists.osgeo.org] *On Behalf Of
*Roxanne Reid-Bennett
*Sent:* 24 November, 2015 12:46 AM
*To:* postgis-users@lists.osgeo.org
<mailto:postgis-users@lists.osgeo.org>
    *Subject:* Re: [postgis-users] Help with SQL query?

On 11/23/2015 12:41 PM, Darrel Maddy wrote:

Dear Pierre,

I was not looking for a total solution and I am grateful for
the suggestion.

Although it may not look like it, I did consult the
documentation and also I have Regina’s book beside me.
 Unfortunately, for me at least, both documents assume some
knowledge of SQL – which I do not have.  I am also trying to
do this simultaneously with a large number of other things
that are new to me.  I do not find the errors reported at all
informative and consider the query I am trying to perform to
be relatively trivial and hence I had hoped the structure of
the query might have been more intuitive.  For others it may be.


FWIW - I don't play with rasters, but this appears to be a pure
SQL thing... add "as rast" like below and try again.


SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF
concentrated > 6 THEN deposition ELSE NULL ENDIF ' )
as rast

FROM mymodel.deposition, mymodel.concentrated

WHERE ST_UpperleftX(mymodel.deposition.rast) =
ST_UpperleftX(mymodel.concentrated.rast) AND

 ST_UpperleftY(mymodel.deposition.rast) =
ST_UpperleftY(mymodel.deposition.rast)

Roxanne




I will try not to bother you again.

Darrel

*From:*postgis-users
[mailto:postgis-users-boun...@lists.osgeo.org] *On Behalf Of
*Pierre Racine
*Sent:* 23 November 2015 20:30
*To:* PostGIS Users Discussion <postgis-users@lists.osgeo.org>
<mailto:postgis-users@lists.osgeo.org>
    *Subject:* Re: [postgis-users] Help with SQL query?

The expression has to stay as it was: 'IF [rast2] > threshold
THEN [rast1] ELSE NULL ENDIF '

Just replace the threshold value as you did.

Do not try to replace the [rast2] and [rast1]. They refer to
the first and second raster pixel values. Read the
ST_Mapalgebra doc…

Don’t expect our suggestions to work blindly. I did not test
this query. I’m not in your context. I expect you read the doc
about all the mentioned functions and adjust for your specific
context. I said “your query should “look like” this”…

Pierre

*From:*postgis-users
[mailto:postgis-users-boun...@lists.osgeo.org] *On Behalf Of
*Darrel Maddy
*Sent:* Monday, November 23, 2015 2:47 PM
*To:* PostGIS Users Discussion
*Subject:* Re: [postgis-users] Help with SQL query?

OK I spoke too soon.

I tried this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum

FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast,
'IF concentrated > 6 THEN deposition ELSE N

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Dear Birgit,

Very many thanks. The query is now running – so it is no longer giving an 
error. I will not know for a few hours whether it is calculating the necessary 
values but it seems progress is being made.

I have a book coming on SQL – I clearly need to do my homework more thoroughly ☺

I will confirm when it stops!

Best wishes

Darrel



From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Birgit Laggner
Sent: 24 November 2015 12:19
To: postgis-users@lists.osgeo.org
Subject: Re: [postgis-users] Help with SQL query?

Hi Darrel,

my PostGIS version is too old for testing, but if I read the documentation 
right, then your expression has to be SQL. And IF THEN ELSE etc. is not SQL as 
far as I know - SQL has CASE WHEN.

So, I would assume, you would need to write your expression like this:

'CASE WHEN [rast2] > 0.0 THEN [rast1] ELSE NULL END'

I am curious if this helps with your error...

Regards,

Birgit

Am 24.11.2015 um 10:09 schrieb Darrel Maddy:
Dear Roxanne,

Many thanks. I did pickup that issue shortly after my last post and it cured 
the first problem.

There is still an SQL issue with the expression however.

I tried ‘IF [rast2]>0.0 THEN [rast1] ELSE NULL ENDIF’

This produces an error at $2. I can remove that problem by doing this
‘IF ([rast2]>0.0) THEN [rast1] ELSE NULL ENDIF’
But then I get the error at THEN

Curiously ‘[rast2]+[rast1]’  works as intended so somehow the conditional is 
being specified incorrectly.

Unfortunately there is nothing in the documentation which helps me with this 
variant.

I will keep trying.

Darrel

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Roxanne Reid-Bennett
Sent: 24 November, 2015 12:46 AM
To: postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

On 11/23/2015 12:41 PM, Darrel Maddy wrote:
Dear Pierre,

I was not looking for a total solution and I am grateful for the suggestion.

Although it may not look like it, I did consult the documentation and also I 
have Regina’s book beside me.  Unfortunately, for me at least, both documents 
assume some knowledge of SQL – which I do not have.  I am also trying to do 
this simultaneously with a large number of other things that are new to me.  I 
do not find the errors reported at all informative and consider the query I am 
trying to perform to be relatively trivial and hence I had hoped the structure 
of the query might have been more intuitive.  For others it may be.

FWIW - I don't play with rasters, but this appears to be a pure SQL thing... 
add "as rast" like below and try again.


SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF concentrated > 6 
THEN deposition ELSE NULL ENDIF ' )
as rast
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast)

Roxanne




I will try not to bother you again.

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 20:30
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org><mailto:postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

The expression has to stay as it was: 'IF [rast2] > threshold THEN [rast1] ELSE 
NULL ENDIF '

Just replace the threshold value as you did.

Do not try to replace the [rast2] and [rast1]. They refer to the first and 
second raster pixel values. Read the ST_Mapalgebra doc…

Don’t expect our suggestions to work blindly. I did not test this query. I’m 
not in your context. I expect you read the doc about all the mentioned 
functions and adjust for your specific context. I said “your query should “look 
like” this”…

Pierre


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 2:47 PM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] Help with SQL query?

OK I spoke too soon.

I tried this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF concentrated 
> 6 THEN deposition ELSE NULL ENDIF ' )
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast) ) foo;

And all I get is rast does not exist.

I’m afraid the penny has not dropped yet ☹

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: 23 November 2015 16:30
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subj

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

OK I spoke too soon.

I tried this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF concentrated 
> 6 THEN deposition ELSE NULL ENDIF ' )
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast) ) foo;

And all I get is rast does not exist.

I'm afraid the penny has not dropped yet :(

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: 23 November 2015 16:30
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

Dear Pierre and Rasmus,

Many thanks for trying to help.

Rasmus: I am aware of GMT but I was looking for a solution in postgis so that I 
can keep all of the data extraction in one place.

Pierre: That is exactly what I was looking for and very many thanks for 
including the explanation.  I am a little overwhelmed with the number of 
functions offered in postgis. It is certainly a remarkable tool.  Watching the 
queries plough through my datasets is a pleasure - albeit the results do not 
always please me :)

Hopefully I can put this to work later tonight.

Best wishes

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 16:17
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,


1)  Create a new raster selecting the right pixels with 
ST_MapAlgebra(raster, raster)

2)  Make sure only intersecting rasters are processed by using their upper 
left corner X and Y coordinates (with ST_UpperLeftX() and UpperLeftY())

3)  Sum the selected pixels with ST_SummaryStats(rast)

All in all a global query should look like this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FFOM (SELECT ST_MapAlgebra(tableA.rast, tableB.rast, 'IF [rast2] > threshold 
THEN [rast1] ELSE NULL ENDIF ' )
FROM tableA, tableB
WHERE ST_UpperleftX(tableA.rast) = ST_UpperleftX(tableB.rast) AND
 ST_UpperleftY(tableA.rast) = 
ST_UpperleftY(tableB.rast) AND
 maybe some other condition here if you get time series 
e.g. tableA.year = tableB.year AND tableA.month = tableB.month) foo

If you have millions of tile you could create indexes on 
ST_UpperleftX(tableA.rast), ST_UpperleftX(tableB.rast), 
ST_UpperleftY(tableA.rast) and ST_UpperleftY(tableB.rast) to make the query 
faster.

You could also just use WHERE ST_Intersects(tableA.rast, tableB.rast) instead...

Pierre

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 7:20 AM
To: PostGIS Users Discussion
Subject: [postgis-users] Help with SQL query?

Dear all,

As you know I am relatively new to postgis and SQL and therefore  I have much 
to learn. However, I am facing a paper deadline and need to do some quick 
analysis of the data I have and I am struggling to figure out how best to 
pursue what I need to do.

I have a significant number of rasters which have double precision values.  
Without going into detail about what the rasters represent, I need to extract 
and sum values from one set of rasters in say table A based upon  values in 
another set of rasters in say table B  where the pixel value in the raster from 
Table B exceeds a threshold. Both tables are the same size (rasters are tiled) 
but I also need to figure out how I make sure the correct rasters are compared. 
 They have filenames like this rastervariable_10.tif, rastervariable_100.tif , 
presumably I need to use a logical expression to strip the numerical value (in 
this case this represents the year) and then order on that basis?

I can do this in QGIS one at a time but that is a little clumsy and rather time 
consuming.

If someone can just point me in the right direction I am sure I can figure out 
the rest for myself.

Apologies once more for asking what is probably a rather trivial question and 
yet again demonstrating my ignorance.

Many thanks

Darrel
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Pierre Racine]

The expression has to stay as it was: 'IF [rast2] > threshold THEN [rast1] ELSE 
NULL ENDIF '

Just replace the threshold value as you did.

Do not try to replace the [rast2] and [rast1]. They refer to the first and 
second raster pixel values. Read the ST_Mapalgebra doc...

Don't expect our suggestions to work blindly. I did not test this query. I'm 
not in your context. I expect you read the doc about all the mentioned 
functions and adjust for your specific context. I said "your query should "look 
like" this"...

Pierre


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 2:47 PM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] Help with SQL query?

OK I spoke too soon.

I tried this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF concentrated 
> 6 THEN deposition ELSE NULL ENDIF ' )
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast) ) foo;

And all I get is rast does not exist.

I'm afraid the penny has not dropped yet :(

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: 23 November 2015 16:30
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Dear Pierre and Rasmus,

Many thanks for trying to help.

Rasmus: I am aware of GMT but I was looking for a solution in postgis so that I 
can keep all of the data extraction in one place.

Pierre: That is exactly what I was looking for and very many thanks for 
including the explanation.  I am a little overwhelmed with the number of 
functions offered in postgis. It is certainly a remarkable tool.  Watching the 
queries plough through my datasets is a pleasure - albeit the results do not 
always please me :)

Hopefully I can put this to work later tonight.

Best wishes

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 16:17
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,


1)  Create a new raster selecting the right pixels with 
ST_MapAlgebra(raster, raster)

2)  Make sure only intersecting rasters are processed by using their upper 
left corner X and Y coordinates (with ST_UpperLeftX() and UpperLeftY())

3)  Sum the selected pixels with ST_SummaryStats(rast)

All in all a global query should look like this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FFOM (SELECT ST_MapAlgebra(tableA.rast, tableB.rast, 'IF [rast2] > threshold 
THEN [rast1] ELSE NULL ENDIF ' )
FROM tableA, tableB
WHERE ST_UpperleftX(tableA.rast) = ST_UpperleftX(tableB.rast) AND
 ST_UpperleftY(tableA.rast) = 
ST_UpperleftY(tableB.rast) AND
 maybe some other condition here if you get time series 
e.g. tableA.year = tableB.year AND tableA.month = tableB.month) foo

If you have millions of tile you could create indexes on 
ST_UpperleftX(tableA.rast), ST_UpperleftX(tableB.rast), 
ST_UpperleftY(tableA.rast) and ST_UpperleftY(tableB.rast) to make the query 
faster.

You could also just use WHERE ST_Intersects(tableA.rast, tableB.rast) instead...

Pierre

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 7:20 AM
To: PostGIS Users Discussion
Subject: [postgis-users] Help with SQL query?

Dear all,

As you know I am relatively new to postgis and SQL and therefore  I have much 
to learn. However, I am facing a paper deadline and need to do some quick 
analysis of the data I have and I am struggling to figure out how best to 
pursue what I need to do.

I have a significant number of rasters which have double precision values.  
Without going into detail about what the rasters represent, I need to extract 
and sum values from one set of rasters in say table A based upon  values in 
another set of rasters in say table B  where the pixel value in the raster from 
Table B exceeds a threshold. Both tables are the same size (rasters are tiled) 
but I also need to figure out how I make sure the correct rasters are compared. 
 They have filenames like this rastervariable_10.tif, rastervariable_100.tif , 
presumably I need to use a logical expression to strip the numerical value (in 
this case this represents the year) and then order on that basis?

I can do this in QGIS one at a time but that is a little clumsy and rather time 
consuming.

If someone can just point me in the right dire

Re: [postgis-users] Help with SQL query?

[Roxanne Reid-Bennett]

On 11/23/2015 12:41 PM, Darrel Maddy wrote:


Dear Pierre,

I was not looking for a total solution and I am grateful for the 
suggestion.


Although it may not look like it, I did consult the documentation and 
also I have Regina’s book beside me.  Unfortunately, for me at least, 
both documents assume some knowledge of SQL – which I do not have.  I 
am also trying to do this simultaneously with a large number of other 
things that are new to me.  I do not find the errors reported at all 
informative and consider the query I am trying to perform to be 
relatively trivial and hence I had hoped the structure of the query 
might have been more intuitive.  For others it may be.




FWIW - I don't play with rasters, but this appears to be a pure SQL 
thing... add "as rast" like below and try again.



SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF 
concentrated > 6 THEN deposition ELSE NULL ENDIF ' )

as rast

FROM mymodel.deposition, mymodel.concentrated

WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND


 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast)



Roxanne



I will try not to bother you again.

Darrel

*From:*postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] 
*On Behalf Of *Pierre Racine

*Sent:* 23 November 2015 20:30
*To:* PostGIS Users Discussion <postgis-users@lists.osgeo.org>
*Subject:* Re: [postgis-users] Help with SQL query?

The expression has to stay as it was: 'IF [rast2] > threshold THEN 
[rast1] ELSE NULL ENDIF '


Just replace the threshold value as you did.

Do not try to replace the [rast2] and [rast1]. They refer to the first 
and second raster pixel values. Read the ST_Mapalgebra doc…


Don’t expect our suggestions to work blindly. I did not test this 
query. I’m not in your context. I expect you read the doc about all 
the mentioned functions and adjust for your specific context. I said 
“your query should “look like” this”…


Pierre

*From:*postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] 
*On Behalf Of *Darrel Maddy

*Sent:* Monday, November 23, 2015 2:47 PM
*To:* PostGIS Users Discussion
*Subject:* Re: [postgis-users] Help with SQL query?

OK I spoke too soon.

I tried this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum

FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF 
concentrated > 6 THEN deposition ELSE NULL ENDIF ' )


FROM mymodel.deposition, mymodel.concentrated

WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND


 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast) ) foo;


And all I get is rast does not exist.

I’m afraid the penny has not dropped yet L

Darrel

*From:*postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] 
*On Behalf Of *Darrel Maddy

*Sent:* 23 November 2015 16:30
*To:* PostGIS Users Discussion <postgis-users@lists.osgeo.org 
<mailto:postgis-users@lists.osgeo.org>>

*Subject:* Re: [postgis-users] Help with SQL query?

Dear Pierre and Rasmus,

Many thanks for trying to help.

Rasmus: I am aware of GMT but I was looking for a solution in postgis 
so that I can keep all of the data extraction in one place.


Pierre: That is exactly what I was looking for and very many thanks 
for including the explanation.  I am a little overwhelmed with the 
number of functions offered in postgis. It is certainly a remarkable 
tool.  Watching the queries plough through my datasets is a pleasure – 
albeit the results do not always please me J


Hopefully I can put this to work later tonight.

Best wishes

Darrel

*From:*postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] 
*On Behalf Of *Pierre Racine

*Sent:* 23 November 2015 16:17
*To:* PostGIS Users Discussion <postgis-users@lists.osgeo.org 
<mailto:postgis-users@lists.osgeo.org>>

*Subject:* Re: [postgis-users] Help with SQL query?

Darrel,

1)Create a new raster selecting the right pixels with 
ST_MapAlgebra(raster, raster)


2)Make sure only intersecting rasters are processed by using their 
upper left corner X and Y coordinates (with ST_UpperLeftX() and 
UpperLeftY())


3)Sum the selected pixels with ST_SummaryStats(rast)

All in all a global query should look like this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum

FFOM (SELECT ST_MapAlgebra(tableA.rast, tableB.rast, 'IF [rast2] > 
threshold THEN [rast1] ELSE NULLENDIF ' )


FROM tableA, tableB

 WHERE ST_UpperleftX(tableA.rast) = ST_UpperleftX(tableB.rast) AND

 ST_UpperleftY(tableA.rast) = 
ST_UpperleftY(tableB.rast) AND


 maybe some other condition here if you get 
time series e.g. tableA.year = tableB.year AND tableA.month = 
tableB.month) foo


If you have millions of tile you could create indexes on 
ST_UpperleftX(tableA.rast), ST_UpperleftX(tableB.rast), 
ST_U

Re: [postgis-users] Help with SQL query?

[Pierre Racine]

Good catch Roxanne! Thanks!

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Roxanne Reid-Bennett
Sent: Monday, November 23, 2015 7:46 PM
To: postgis-users@lists.osgeo.org
Subject: Re: [postgis-users] Help with SQL query?

On 11/23/2015 12:41 PM, Darrel Maddy wrote:
Dear Pierre,

I was not looking for a total solution and I am grateful for the suggestion.

Although it may not look like it, I did consult the documentation and also I 
have Regina’s book beside me.  Unfortunately, for me at least, both documents 
assume some knowledge of SQL – which I do not have.  I am also trying to do 
this simultaneously with a large number of other things that are new to me.  I 
do not find the errors reported at all informative and consider the query I am 
trying to perform to be relatively trivial and hence I had hoped the structure 
of the query might have been more intuitive.  For others it may be.

FWIW - I don't play with rasters, but this appears to be a pure SQL thing... 
add "as rast" like below and try again.


SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF concentrated > 6 
THEN deposition ELSE NULL ENDIF ' )
as rast
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast)

Roxanne



I will try not to bother you again.

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 20:30
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org><mailto:postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

The expression has to stay as it was: 'IF [rast2] > threshold THEN [rast1] ELSE 
NULL ENDIF '

Just replace the threshold value as you did.

Do not try to replace the [rast2] and [rast1]. They refer to the first and 
second raster pixel values. Read the ST_Mapalgebra doc…

Don’t expect our suggestions to work blindly. I did not test this query. I’m 
not in your context. I expect you read the doc about all the mentioned 
functions and adjust for your specific context. I said “your query should “look 
like” this”…

Pierre


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 2:47 PM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] Help with SQL query?

OK I spoke too soon.

I tried this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF concentrated 
> 6 THEN deposition ELSE NULL ENDIF ' )
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast) ) foo;

And all I get is rast does not exist.

I’m afraid the penny has not dropped yet ☹

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: 23 November 2015 16:30
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Dear Pierre and Rasmus,

Many thanks for trying to help.

Rasmus: I am aware of GMT but I was looking for a solution in postgis so that I 
can keep all of the data extraction in one place.

Pierre: That is exactly what I was looking for and very many thanks for 
including the explanation.  I am a little overwhelmed with the number of 
functions offered in postgis. It is certainly a remarkable tool.  Watching the 
queries plough through my datasets is a pleasure – albeit the results do not 
always please me ☺

Hopefully I can put this to work later tonight.

Best wishes

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 16:17
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,


1)  Create a new raster selecting the right pixels with 
ST_MapAlgebra(raster, raster)

2)  Make sure only intersecting rasters are processed by using their upper 
left corner X and Y coordinates (with ST_UpperLeftX() and UpperLeftY())

3)  Sum the selected pixels with ST_SummaryStats(rast)

All in all a global query should look like this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FFOM (SELECT ST_MapAlgebra(tableA.rast, tableB.rast, 'IF [rast2] > threshold 
THEN [rast1] ELSE NULL ENDIF ' )
FROM tableA, tableB
WHERE ST_UpperleftX(tableA.rast) = ST_UpperleftX(tableB.rast) AND
 ST_UpperleftY(tableA.rast) = 
ST_Upper

Re: [postgis-users] Help with SQL query?

[Rasmus Aveskogh]

Regarding filename matching etc. I would just wrap everything in a scripting 
language of your choice (bash, python, ruby etc.) considering all the GMT 
programs, including grdmath, are command line tools.

-ra

On 23 Nov 2015, at 16:27, Rasmus Aveskogh  wrote:

> 
> Though this is probably possible in PostGIS nowadays I would personally solve 
> such a use case with GMT (http://gmt.soest.hawaii.edu/), for example with 
> grdmath (http://gmt.soest.hawaii.edu/doc/5.2.1/grdmath.html) which lets you 
> perform rather complex computations on multiple grids (for examples rasters) 
> using reverse polish notation.
> 
> -ra
> 
> On 23 Nov 2015, at 13:20, Darrel Maddy  wrote:
> 
>> Dear all,
>>  
>> As you know I am relatively new to postgis and SQL and therefore  I have 
>> much to learn. However, I am facing a paper deadline and need to do some 
>> quick analysis of the data I have and I am struggling to figure out how best 
>> to pursue what I need to do.
>>  
>> I have a significant number of rasters which have double precision values.  
>> Without going into detail about what the rasters represent, I need to 
>> extract and sum values from one set of rasters in say table A based upon  
>> values in another set of rasters in say table B  where the pixel value in 
>> the raster from Table B exceeds a threshold. Both tables are the same size 
>> (rasters are tiled) but I also need to figure out how I make sure the 
>> correct rasters are compared.  They have filenames like this 
>> rastervariable_10.tif, rastervariable_100.tif , presumably I need to use a 
>> logical expression to strip the numerical value (in this case this 
>> represents the year) and then order on that basis?
>>  
>> I can do this in QGIS one at a time but that is a little clumsy and rather 
>> time consuming.
>>  
>> If someone can just point me in the right direction I am sure I can figure 
>> out the rest for myself.
>>  
>> Apologies once more for asking what is probably a rather trivial question 
>> and yet again demonstrating my ignorance.
>>  
>> Many thanks
>>  
>> Darrel
>> ___
>> postgis-users mailing list
>> postgis-users@lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/postgis-users
> 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Rasmus Aveskogh]

Though this is probably possible in PostGIS nowadays I would personally solve 
such a use case with GMT (http://gmt.soest.hawaii.edu/), for example with 
grdmath (http://gmt.soest.hawaii.edu/doc/5.2.1/grdmath.html) which lets you 
perform rather complex computations on multiple grids (for examples rasters) 
using reverse polish notation.

-ra

On 23 Nov 2015, at 13:20, Darrel Maddy  wrote:

> Dear all,
>  
> As you know I am relatively new to postgis and SQL and therefore  I have much 
> to learn. However, I am facing a paper deadline and need to do some quick 
> analysis of the data I have and I am struggling to figure out how best to 
> pursue what I need to do.
>  
> I have a significant number of rasters which have double precision values.  
> Without going into detail about what the rasters represent, I need to extract 
> and sum values from one set of rasters in say table A based upon  values in 
> another set of rasters in say table B  where the pixel value in the raster 
> from Table B exceeds a threshold. Both tables are the same size (rasters are 
> tiled) but I also need to figure out how I make sure the correct rasters are 
> compared.  They have filenames like this rastervariable_10.tif, 
> rastervariable_100.tif , presumably I need to use a logical expression to 
> strip the numerical value (in this case this represents the year) and then 
> order on that basis?
>  
> I can do this in QGIS one at a time but that is a little clumsy and rather 
> time consuming.
>  
> If someone can just point me in the right direction I am sure I can figure 
> out the rest for myself.
>  
> Apologies once more for asking what is probably a rather trivial question and 
> yet again demonstrating my ignorance.
>  
> Many thanks
>  
> Darrel
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Pierre Racine]

Darrel,


1)  Create a new raster selecting the right pixels with 
ST_MapAlgebra(raster, raster)

2)  Make sure only intersecting rasters are processed by using their upper 
left corner X and Y coordinates (with ST_UpperLeftX() and UpperLeftY())

3)  Sum the selected pixels with ST_SummaryStats(rast)

All in all a global query should look like this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FFOM (SELECT ST_MapAlgebra(tableA.rast, tableB.rast, 'IF [rast2] > threshold 
THEN [rast1] ELSE NULL ENDIF ' )
FROM tableA, tableB
WHERE ST_UpperleftX(tableA.rast) = ST_UpperleftX(tableB.rast) AND
 ST_UpperleftY(tableA.rast) = 
ST_UpperleftY(tableB.rast) AND
 maybe some other condition here if you get time series 
e.g. tableA.year = tableB.year AND tableA.month = tableB.month) foo

If you have millions of tile you could create indexes on 
ST_UpperleftX(tableA.rast), ST_UpperleftX(tableB.rast), 
ST_UpperleftY(tableA.rast) and ST_UpperleftY(tableB.rast) to make the query 
faster.

You could also just use WHERE ST_Intersects(tableA.rast, tableB.rast) instead...

Pierre

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 7:20 AM
To: PostGIS Users Discussion
Subject: [postgis-users] Help with SQL query?

Dear all,

As you know I am relatively new to postgis and SQL and therefore  I have much 
to learn. However, I am facing a paper deadline and need to do some quick 
analysis of the data I have and I am struggling to figure out how best to 
pursue what I need to do.

I have a significant number of rasters which have double precision values.  
Without going into detail about what the rasters represent, I need to extract 
and sum values from one set of rasters in say table A based upon  values in 
another set of rasters in say table B  where the pixel value in the raster from 
Table B exceeds a threshold. Both tables are the same size (rasters are tiled) 
but I also need to figure out how I make sure the correct rasters are compared. 
 They have filenames like this rastervariable_10.tif, rastervariable_100.tif , 
presumably I need to use a logical expression to strip the numerical value (in 
this case this represents the year) and then order on that basis?

I can do this in QGIS one at a time but that is a little clumsy and rather time 
consuming.

If someone can just point me in the right direction I am sure I can figure out 
the rest for myself.

Apologies once more for asking what is probably a rather trivial question and 
yet again demonstrating my ignorance.

Many thanks

Darrel
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Dear Pierre and Rasmus,

Many thanks for trying to help.

Rasmus: I am aware of GMT but I was looking for a solution in postgis so that I 
can keep all of the data extraction in one place.

Pierre: That is exactly what I was looking for and very many thanks for 
including the explanation.  I am a little overwhelmed with the number of 
functions offered in postgis. It is certainly a remarkable tool.  Watching the 
queries plough through my datasets is a pleasure - albeit the results do not 
always please me :)

Hopefully I can put this to work later tonight.

Best wishes

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 16:17
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,


1)  Create a new raster selecting the right pixels with 
ST_MapAlgebra(raster, raster)

2)  Make sure only intersecting rasters are processed by using their upper 
left corner X and Y coordinates (with ST_UpperLeftX() and UpperLeftY())

3)  Sum the selected pixels with ST_SummaryStats(rast)

All in all a global query should look like this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FFOM (SELECT ST_MapAlgebra(tableA.rast, tableB.rast, 'IF [rast2] > threshold 
THEN [rast1] ELSE NULL ENDIF ' )
FROM tableA, tableB
WHERE ST_UpperleftX(tableA.rast) = ST_UpperleftX(tableB.rast) AND
 ST_UpperleftY(tableA.rast) = 
ST_UpperleftY(tableB.rast) AND
 maybe some other condition here if you get time series 
e.g. tableA.year = tableB.year AND tableA.month = tableB.month) foo

If you have millions of tile you could create indexes on 
ST_UpperleftX(tableA.rast), ST_UpperleftX(tableB.rast), 
ST_UpperleftY(tableA.rast) and ST_UpperleftY(tableB.rast) to make the query 
faster.

You could also just use WHERE ST_Intersects(tableA.rast, tableB.rast) instead...

Pierre

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 7:20 AM
To: PostGIS Users Discussion
Subject: [postgis-users] Help with SQL query?

Dear all,

As you know I am relatively new to postgis and SQL and therefore  I have much 
to learn. However, I am facing a paper deadline and need to do some quick 
analysis of the data I have and I am struggling to figure out how best to 
pursue what I need to do.

I have a significant number of rasters which have double precision values.  
Without going into detail about what the rasters represent, I need to extract 
and sum values from one set of rasters in say table A based upon  values in 
another set of rasters in say table B  where the pixel value in the raster from 
Table B exceeds a threshold. Both tables are the same size (rasters are tiled) 
but I also need to figure out how I make sure the correct rasters are compared. 
 They have filenames like this rastervariable_10.tif, rastervariable_100.tif , 
presumably I need to use a logical expression to strip the numerical value (in 
this case this represents the year) and then order on that basis?

I can do this in QGIS one at a time but that is a little clumsy and rather time 
consuming.

If someone can just point me in the right direction I am sure I can figure out 
the rest for myself.

Apologies once more for asking what is probably a rather trivial question and 
yet again demonstrating my ignorance.

Many thanks

Darrel
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Darrel Maddy]

Dear Pierre,

I was not looking for a total solution and I am grateful for the suggestion.

Although it may not look like it, I did consult the documentation and also I 
have Regina's book beside me.  Unfortunately, for me at least, both documents 
assume some knowledge of SQL - which I do not have.  I am also trying to do 
this simultaneously with a large number of other things that are new to me.  I 
do not find the errors reported at all informative and consider the query I am 
trying to perform to be relatively trivial and hence I had hoped the structure 
of the query might have been more intuitive.  For others it may be.

I will try not to bother you again.

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 20:30
To: PostGIS Users Discussion <postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

The expression has to stay as it was: 'IF [rast2] > threshold THEN [rast1] ELSE 
NULL ENDIF '

Just replace the threshold value as you did.

Do not try to replace the [rast2] and [rast1]. They refer to the first and 
second raster pixel values. Read the ST_Mapalgebra doc...

Don't expect our suggestions to work blindly. I did not test this query. I'm 
not in your context. I expect you read the doc about all the mentioned 
functions and adjust for your specific context. I said "your query should "look 
like" this"...

Pierre


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 2:47 PM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] Help with SQL query?

OK I spoke too soon.

I tried this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF concentrated 
> 6 THEN deposition ELSE NULL ENDIF ' )
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast) ) foo;

And all I get is rast does not exist.

I'm afraid the penny has not dropped yet :(

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: 23 November 2015 16:30
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Dear Pierre and Rasmus,

Many thanks for trying to help.

Rasmus: I am aware of GMT but I was looking for a solution in postgis so that I 
can keep all of the data extraction in one place.

Pierre: That is exactly what I was looking for and very many thanks for 
including the explanation.  I am a little overwhelmed with the number of 
functions offered in postgis. It is certainly a remarkable tool.  Watching the 
queries plough through my datasets is a pleasure - albeit the results do not 
always please me :)

Hopefully I can put this to work later tonight.

Best wishes

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 16:17
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,


1)  Create a new raster selecting the right pixels with 
ST_MapAlgebra(raster, raster)

2)  Make sure only intersecting rasters are processed by using their upper 
left corner X and Y coordinates (with ST_UpperLeftX() and UpperLeftY())

3)  Sum the selected pixels with ST_SummaryStats(rast)

All in all a global query should look like this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FFOM (SELECT ST_MapAlgebra(tableA.rast, tableB.rast, 'IF [rast2] > threshold 
THEN [rast1] ELSE NULL ENDIF ' )
FROM tableA, tableB
WHERE ST_UpperleftX(tableA.rast) = ST_UpperleftX(tableB.rast) AND
 ST_UpperleftY(tableA.rast) = 
ST_UpperleftY(tableB.rast) AND
 maybe some other condition here if you get time series 
e.g. tableA.year = tableB.year AND tableA.month = tableB.month) foo

If you have millions of tile you could create indexes on 
ST_UpperleftX(tableA.rast), ST_UpperleftX(tableB.rast), 
ST_UpperleftY(tableA.rast) and ST_UpperleftY(tableB.rast) to make the query 
faster.

You could also just use WHERE ST_Intersects(tableA.rast, tableB.rast) instead...

Pierre

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 7:20 AM
To: PostGIS Users Discussion
Subject: [postgis-users] Help with SQL query?

Dear all,

As you know I am relatively new to postgis and SQL and therefore  I have much 
to learn. However, I am facing a paper deadline and need to do some quick 
analysis of