[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-bounces at lists.osgeo.org] On Behalf Of Darrel Maddy
Sent: 28 November 2015 22:51
To: PostGIS Users Discussion <postgis-users at 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-bounces at lists.osgeo.org] On Behalf Of Paragon Corporation
Sent: 28 November 2015 05:04
To: 'PostGIS Users Discussion' <postgis-users at lists.osgeo.org<mailto:postgis-users at 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-10000]: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 10000 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-10000 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-bounces at lists.osgeo.org] On Behalf Of Paragon Corporation
Sent: Friday, November 27, 2015 10:04 PM
To: 'PostGIS Users Discussion' <postgis-users at lists.osgeo.org<mailto:postgis-users at 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-bounces at lists.osgeo.org] On Behalf Of Darrel Maddy
Sent: Friday, November 27, 2015 4:09 PM
To: PostGIS Users Discussion <postgis-users at lists.osgeo.org<mailto:postgis-users at 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-bounces at lists.osgeo.org] On Behalf Of Paragon Corporation
Sent: 27 November 2015 20:02
To: 'PostGIS Users Discussion' <postgis-users at lists.osgeo.org<mailto:postgis-users at 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20151128/0e6f760f/attachment.html>