Re: [postgis-users] Incorrect Result for ST_CountAgg with nodata pixels?

[Michael Treglia]

Hi All,

Just to follow-up - looks like I've figured out the issue:
In the count of total number of pixels per building, it seems I was 
calculating the total number of pixels per tile that intersected with the 
polygons, not the total number of pixels that actually intersected the 
polygon. For now, I used a reclassification step, which gets around this 
issue, but I suspect an ST_Intersection approach would also work. 

Here's the working sql:

--Begin Code
SELECT *,
 count_undevel/(count_total)::float as PropVeg --Calculates proportion of 
pixels within each building that are class 1 or 2 based on count_undevel 
and count_total
INTO resultsvectorlayers.greenroofs_pctgrn
FROM(
SELECT 
p.*,
--Count number of pixels per building with values 1 or 2
ST_CountAgg(
ST_Reclass(
ST_Clip(rast, p.geom_2263),
1,'[1-2]:1, [3-5]:0', '8BUI',0), 1,TRUE) --reclassifies classes 1 and 2 as 
1, reclassifies values 3-5 as 0, and designates 0 as nodata
as count_undevel,
as count_undevel,
--Count total number of pixels
ST_CountAgg(
ST_Clip(
ST_Reclass(rast,1,'[0-5]:1','8BUI',1), --reclassify all pixels to 
non-nodata; when this is clipped, all pixels within polygons are kept as 1, 
while pixels outside of buildings are set to nodata
1, p.geom_2263,NULL,TRUE),
1,TRUE) --exclude nodata from countAgg
as count_total
FROM staging.naip_ndvi_classified_2013 r, staging.nycbldgs p
WHERE ST_INTERSECTS(r.rast, p.geom_2263) --and rid = 2
group by p.geom_2263, p.gid
) AS FOO;
--End Code

On Tuesday, July 12, 2016 at 4:56:11 PM UTC-4, Michael Treglia wrote:
>
> Hi All,
>
> I'm working to do some overlay analyses between raster and vector layers, 
> and having an issue getting an accurate count of pixels per polygon.
>
> My query is below - basically my raster layer has values ranging 1-5, with 
> some nodata pixels (defined as such, with values of 0). I want to 
> calculate, for each polygon, the proportion of pixels, of the total 
> (including nodata), with values of 1 and 2. 
>
> The first ST_CountAgg statement (yielding 'count_undevel') seems to be 
> working as expected, and reports an accurate number. However, the second 
> one, to create the 'count_total' field consistently gives an over-count.  
>
> Any guess as to what might be going on/how to correct? I'll gladly share 
> the properties of the data if that would be useful.
>
> Thanks!
> Mike
>
>
>
> SELECT *,
>  count_undevel/(count_total)::float as PropVeg --Calculates proportion of 
> pixels within each building that are class 1 or 2 based on count_undevel 
> and count_total
> INTO resultsvectorlayers.pctgrn
> FROM(
> SELECT 
> p.*,
> --Count number of pixels per building with values 1 or 2
> ST_CountAgg(
> ST_Reclass(
> ST_Clip(rast, p.geom_2263),
> 1,'[1-2]:1, [3-5]:0', '8BUI',0), 1,TRUE) --reclassifies classes 1 and 2 as 
> 1, reclassifies values 3-5 as 0, and designates 0 as nodata
> as count_undevel,
> --Count total number of pixels
> ST_CountAgg(
> ST_Clip(rast, p.geom_2263),
> 1,FALSE)
> as count_total
> FROM staging.naip_ndvi_classified_2013 r, staging.nycbldgs p
> WHERE ST_INTERSECTS(r.rast, p.geom_2263) --and rid = 2
> group by p.geom_2263, p.gid
> ) AS FOO;
>
>
>
>___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] number all houses on branch roads first

[Rémi Cura]

Would be much easier to send your graph to python networkx

If you have plpython, this is immediate

(line 368)

I found it painful  to work on graph in pgsql (or pgrouting, besides the
excellent already written functions),

much easier with dedicated tools.


If you really really want to use pgsql,
you can use pure SQL with a  recursive CTE,
as in postgis

topology.


I find it hard to test and debug though

Cheers,
Rémi-C

2016-07-13 4:17 GMT+02:00 Stephen Woodbridge :

> Opened an enhancement request to pgrouting:
> https://github.com/pgRouting/pgrouting/issues/625
>
>
> On 7/12/2016 9:56 PM, Stephen Woodbridge wrote:
>
>> Sounds like you want a "graph" based on the road network topology like
>> what is used in pgRouting. Then given the graph do a depth first
>> traversal of the topology, labeling the edges as you go.
>>
>> If you build the graph topology (not the postgis topology), the you can
>> probably write a recursive query to do the traversal and labeling. I
>> think I would probably write a pgsql function(s) to traverse the graph
>> and label the edges.
>>
>> Unfortunately, pgRouting does not have depth first search function built
>> into it, hence the need to code one in pgsql.
>>
>> -Steve
>>
>> On 7/12/2016 5:16 PM, Dan Jacobson wrote:
>>
>>> I want to go down the (all unnamed) roads in my future mountain
>>> community assigning house numbers every 50 meters.
>>>
>>>50-100   170-200
>>> / /
>>> 1  / 107 147 /  213
>>> +-+---+---+-++--240--main-road--
>>> 2 48   \   168
>>> \
>>>   120-140
>>>
>>> I stay on my main road, but whenever encountering a fork, first go
>>> down it.
>>> "Depth first pre-order ordered labeled rooted binary tree traversal
>>> but with central path"?
>>>
>>> http://math.stackexchange.com/questions/1856814/binary-tree-traversal-with-fixed-final-node
>>>
>>>
>>> OK for PostGIS, given a few linestrings, I suppose I first (somehow)
>>> connect them to form a network, then ride my virtual car down it. And
>>> whenever my odometer reaches another 50m, make a mark on the centerline.
>>> (Assume a strict binary tree (mountain roads with no 4-way junctions))
>>> At each road junction I first choose a side road, and if already on a
>>> side road, first choose the left road, before choosing the right road.
>>> When backtracking turn off my odometer, until finally back on to my main
>>> road.
>>>
>>> Holy smokes, sounds tough. Can I do this with PostGIS or should I go
>>> back to GRASS or what?
>>>
>>> (I am thinking instead of using ST_OffsetCurve (previous project, thanks
>>> Sandro) to put odd on the left even on the right, down the centerline
>>> I'll just put odd numbers at 25, 75m, padded "1 ", and even at 50, 100m
>>> reverse padded " 2".) These points and their labels are what I want
>>> for output.
>>> ___
>>> postgis-users mailing list
>>> postgis-users@lists.osgeo.org
>>> http://lists.osgeo.org/mailman/listinfo/postgis-users
>>>
>>>
>>
>> ---
>> This email has been checked for viruses by Avast antivirus software.
>> https://www.avast.com/antivirus
>>
>> ___
>> postgis-users mailing list
>> postgis-users@lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/postgis-users
>>
>
>
> ---
> This email has been checked for viruses by Avast antivirus software.
> https://www.avast.com/antivirus
>
> ___
> 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] Question on topology

[Sandro Santilli]

On Wed, Jul 13, 2016 at 01:05:41PM +0200, Neumann, Andreas wrote:

> Yes - I will try that next - loading it in chunks. But wouldn't it miss
> out on some of the neighbourpolygons then if I use such subsets based on
> pkey? 

The final goal is importing all of your geometries, so eventually
neighbours would be added soon or late.

If you prefer loading in order you could use use
ORDER BY geom OFFSET $1 LIMIT $2, which would insert
inputs from left to right ...

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

Re: [postgis-users] Question on topology

[Neumann, Andreas]

Hi Strk, 

Yes - I will try that next - loading it in chunks. But wouldn't it miss
out on some of the neighbourpolygons then if I use such subsets based on
pkey? 

Thanks,
Andreas 

On 2016-07-13 12:55, Sandro Santilli wrote:

> On Tue, Jul 12, 2016 at 01:15:41PM +0200, Neumann, Andreas wrote: 
> 
>> Hi Sandro, 
>> 
>> I let this function run during the night and after 7.5h I got an error
>> message:
> 
> Eh, that's why I suggested loading the topology in batches.
> 
> Use constructs like:
> 
> SELECT TopoGeo_addPolygon(...)
> FROM input_table
> WHERE gid >= $1 AND gid < $2
> 
> With $1..$2 in a range that makes it run within a few minutes.
> 
> Or, alternatively, code the loop in plpgsql by intercepting
> exceptions and skipping the "offending" input to analyze later
> (but I like the multi-transaction approach better as it allows
> you to see the topology in QGIS while it's being loaded).
> 
>> ERROR: Corrupted topology: adjacent edges 159972 and -159958 bind
>> different face (0 and 78697)
> 
> It means that somehow the ST_CreateTopoGeo function in a given
> stage created a corrupted topology. This is usually due to some
> robustness issue. Snapping may help a little with this.
> 
> Once again, note that if this was being done in a loop you might
> have stopped and looked at the issue with QGIS, to continue after
> the topology validity was fixed.
> 
>> I have to say that I did not tune Postgis memory-wise - maybe I should
>> do that next. But the above error message probably doesn't indicate a
>> memory problem, but a problem with the data.
> 
> Correct, it's not a problem with memory.
> 
>> My data: 2582 polygons with a total of 167176 vertices.
> 
> Try loading it in chunks of 300 polygons
> 
> --strk;
> 
> ()   Free GIS & Flash consultant/developer
> /\   https://strk.kbt.io/services.html

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