[postgis-users] Declarative partitioning by hash on geom column

[Birgit Laggner]

Dear list members, 

I have a table with 35 billion rows containing point geometries (SRID 31467) 
and use it for intersections with multiple polygon datasets. To improve 
performance I would like to use declarative partitioning to split the point 
table into sensible table partitions. Because of the intended use with 
ST_Intersects, I thought it would somehow be best to use geometric location as 
the basis for partitioning. List and range partitioning are not really 
applicable to the geom column, so that leaves hash partitioning. But I get the 
feeling that hash partioning won't work for the geometry data type either. 

All posts I found on the internet about partitioning PostGIS tables were mostly 
outdated and used inheritance partioning with check constraints. 

Has anyone used declarative partitioning on PostGIS geometry tables or can 
otherwise tell me why this might be impossible without using an additional 
column? 

Help would be very much appreciated :-) 

Thanks a lot and kind regards 
Birgit 

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

Re: [postgis-users] Location of highest point in each region?

[Birgit Laggner]

Hi Jesper, 

just from the manual: If you have the value, could you not use ST_PixelOfValue 
to extract the location? 

from the PostGIS-manual: 



ST_PixelOfValue — Get the columnx, rowy coordinates of the pixel whose value 
equals the search value. 



Regards, 
Birgit 



Von: "Jesper Kempe"  
An: postgis-users@lists.osgeo.org 
Gesendet: Montag, 19. April 2021 08:34:42 
Betreff: [postgis-users] Location of highest point in each region? 

Hi 
I have loaded a DEM dataset into a Postgis raster table. And a shapefile with 
polygons for a number of regions into another table. 

For a couple of days now I have been trying to extract the height and location 
of the highest point in each of the regions. But I cannot find a way to do it. 

I am able to find the highest value in each region using ST_SummaryStats() but 
not the location (coordinates) of that value. 

Can anyone point me in the right direction how I should do this? 
It must be possible to do this in Postgis right? 

Thanks 
Jesper 

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

Re: [postgis-users] Stuck with PostgreSQL/PostGIS Upgrade

[Birgit Laggner]

Dear Paul, 

thank you very much for your helpful response! We will try the symlink and I 
will read your blog to learn more about the upgrade process. And yes, sorry for 
the limited information, but it was already quite difficult to get these 
snippets from our admin... :-) 

With thankful regards, 
Birgit 



Von: "Paul Ramsey"  
An: "PostGIS Users Discussion"  
Gesendet: Dienstag, 13. August 2019 17:46:42 
Betreff: Re: [postgis-users] Stuck with PostgreSQL/PostGIS Upgrade 

Some more background information about the magic gears turning behind the 
upgrade process: 

[ http://blog.cleverelephant.ca/2016/08/postgis-upgrade.html | 
http://blog.cleverelephant.ca/2016/08/postgis-upgrade.html ] 

P 




On Aug 13, 2019, at 8:40 AM, Paul Ramsey < [ mailto:pram...@cleverelephant.ca | 
pram...@cleverelephant.ca ] > wrote: 

B, 
You haven’t included quite enough information to make a full diagnosis (like 
what command generates the error message) but just guessing, your admin used 
apt-get to update both postgres and postgis packages and then ran pg_upgrade on 
the cluster. Since the old postgis version is gone at this point (it’s been 
removed during the update) when the system catalog tries to make sense of 
postgis objects and references to old library, you get this error. 
A hack to work around the issue is to find the location of postgis-2.5.so (the 
new library copy from the upgrade) and then create a symlink from the 
postgis-2.3.so to that file, so that the system will end up finding a library 
to load. The new and old libraries have mostly the same symbols, so this trick 
works fine for the purposes of upgrade. 
You may run into yet further problems, but this is a quick way to try and move 
forward. Once you have the database up and running, do not forget to run ALTER 
'EXTENSION postgis UDPATE’ on all your postgis databases. This re-installs all 
the SQL functions, pointing to the new library file. Once that’s done you can 
remove the symlink. 
ATB, 
P 

PS, here’s a link to a blog post about working through this same issue, but on 
a Centos system. So it’s not command-perfect, but provides general explanations 
of problems that an admin might use to navigate the Ubuntu specific issues. [ 
https://info.crunchydata.com/blog/upgrading-postgis-on-centos-7 | 
https://info.crunchydata.com/blog/upgrading-postgis-on-centos-7 ] 


BQ_BEGIN
On Aug 13, 2019, at 7:56 AM, Birgit Laggner < [ 
mailto:birgit.lagg...@thuenen.de | birgit.lagg...@thuenen.de ] > wrote: 

Dear list, 

on our test server, we currently try to test how to upgrade our 
PostgreSQL/PostGIS database. Sadly, we are stuck somewhere and would be very 
glad if anybody could help. 

Our test server is running on an Ubuntu 16.04 system in a virtual machine. We 
are coming from PostgreSQL 9.5.17 with PostGIS 2.3.3 and want to upgrade to 
PostgreSQL 11.2 and PostGIS 2.5. 

Here the outputs of the version queries: 
SELECT version(); 
"PostgreSQL 9.5.17 on x86_64-pc-linux-gnu, compiled by gcc (Ubuntu 
5.4.0-6ubuntu1~16.04.10) 5.4.0 20160609, 64-bit" 

SELECT postgis_full_version(); 
"POSTGIS="2.3.3 r15473" GEOS="3.5.1-CAPI-1.9.1 r4246" SFCGAL="1.2.2" PROJ="Rel. 
4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 2015/09/16" 
LIBXML="2.9.3" LIBJSON="0.11.99" TOPOLOGY RASTER" 

Our admin tried to upgrade PostGIS first and swears he did as told in the 
PostGIS manual. But apparently, the upgrade was not fully successful, since 
PostgreSQL still seems to look for PostGIS 2.3. 

His error message while trying to upgrade PostgreSQL is: 
postgres@gis ERROR: could not access file "$libdir/postgis-2.3": No such file 
or directory. 

Has anybody any idea what we are doing wrong and what we should be doing 
instead? 

With hopeful regards, 
Birgit 


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




BQ_END



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

[postgis-users] Stuck with PostgreSQL/PostGIS Upgrade

[Birgit Laggner]

Dear list, 

on our test server, we currently try to test how to upgrade our 
PostgreSQL/PostGIS database. Sadly, we are stuck somewhere and would be very 
glad if anybody could help. 

Our test server is running on an Ubuntu 16.04 system in a virtual machine. We 
are coming from PostgreSQL 9.5.17 with PostGIS 2.3.3 and want to upgrade to 
PostgreSQL 11.2 and PostGIS 2.5. 

Here the outputs of the version queries: 
SELECT version(); 
"PostgreSQL 9.5.17 on x86_64-pc-linux-gnu, compiled by gcc (Ubuntu 
5.4.0-6ubuntu1~ 16.04.10 ) [ callto:5.4.0 20160609, 64 | 5.4.0 20160609, 64 ] 
-bit" 

SELECT postgis_full_version(); 
"POSTGIS="2.3.3 r15473" GEOS="3.5.1-CAPI-1.9.1 r4246" SFCGAL="1.2.2" PROJ="Rel. 
4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 2015/09/16" 
LIBXML="2.9.3" LIBJSON="0.11.99" TOPOLOGY RASTER" 

Our admin tried to upgrade PostGIS first and swears he did as told in the 
PostGIS manual. But apparently, the upgrade was not fully successful, since 
PostgreSQL still seems to look for PostGIS 2.3. 

His error message while trying to upgrade PostgreSQL is: 
postgres@gis ERROR: could not access file "$libdir/postgis-2.3": No such file 
or directory. 

Has anybody any idea what we are doing wrong and what we should be doing 
instead? 

With hopeful regards, 
Birgit 


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

Re: [postgis-users] Geography/geometry

[Birgit Laggner]

Hi Simon, 

I assume, you are trying to get this all done in one query, right? Did you try 
to throw ST_MakeValid into the mix? Like: 
ST_Transform(ST_Union(ST_MakeValid(ST_Transform(a.geog4326::geometry, 3857))), 
4326)::geography 

Regards, 
Birgit 




Von: "Simon Greener"  
An: "PostGIS Users Discussion"  
Gesendet: Montag, 21. Januar 2019 06:23:10 
Betreff: [postgis-users] Geography/geometry 

Folks, 

I'm in a situation where I need to ST_Union or ST_Collect some osm_county 
Polygon (not MultiPolygon) data for Ireland. 

Now, because ST_Union or ST_Collect do not support geography, I cast to 
geometry before calling. 

select min(a.osm_id) as osm_id, 
a.name, 
count(*) as parts, 
ST_Union(a.geog4326::geometry)::geography as geog4326 -- or ST_Collect 
from data.osm_county as a 
group by a.name; 

Whence I get this: 

ERROR: lwgeom_area_spher(oid) returned area < 0.0 

Investigating I get results like this: 

select distinct st_isvalidreason(a.geog4326::geometry) from data.osm_county as 
a; 

"Hole lies outside shell[-10.2589459 53.9746452]" 
etc 

I guess this is expected because geodetic lines in the source geography are 
being treated as straight in the cast'd geometry. 

If I use ST_Transform to project a 4326 poly to a 3857 and then call the 
ST_Union aggregate, or identify a single geography that has the invalidity and 
execute a self-union, I get the following in both situations. 

ERROR: GEOSUnaryUnion: TopologyException: Input geom 0 is invalid: Hole lies 
outside shell at or near point -1148162.9982628345 7095296.1166736316 at 
-1148162.9982628345 7095296.1166736316 

I can't for the life of me work out how to complete the aggregated ST_Union on 
the 4326 geography data. 

Anyone point out what I am doing wrong or give me a pointer to what I can do to 
achieve the aggregated union? 
Regards 
Simon 

 
Spatial Advice & Solutions Architecture 
Database Spatial Stored Procedure Designer 
Oracle Spatial, SQL Server, PostGIS, MySQL, ArcSDE FME 
Awarded "2011 Oracle Spatial Excellence Award for Education and Research" 
A: 39 Cliff View Drive, Allens Rivulet, 7150, Tas, Aust 
W: www.spdba.com.au 
E: si...@spdba.com.au 
V: +61 362 396 397 
M: +61 418 396 391 
GITC Supplier: T1005 
Skype: sggreener 
Long: 147.20515 (147° 12' 18" E) 
Lat: -43.01530 (43° 00' 55" S) 
GeoHash: r22em9r98wg 
NAC:W80CK 7SWP3 

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

Re: [postgis-users] PostGIS case usages

[Birgit Laggner]

Hi Regina,

I work at a Federal Research Institute and we use PostGIS for various purposes. 
A very common use case is joining different spatial datasets (vector and 
sometimes raster) with ST_Intersects, e.g. to find out how agricultural land 
use is related to nature conservation areas or other protection areas, or to 
determine the land use changes over time to evaluate the effectiveness of 
certain political measures.

As far as I know, most PostGIS users at our institute have PostgreSQL databases 
on Linux machines.

Does that help?

Regards and success with the new book!

Birgit

- Ursprüngliche Mail -
Von: "Regina Obe" 
An: postgis-users@lists.osgeo.org
CC: "PostGIS Development Discussion" 
Gesendet: Freitag, 26. Oktober 2018 19:05:52
Betreff: [postgis-users] PostGIS case usages

Hey all.  So we've been in talks with our editor about having a 3rd Edition
of PostGIS hopefully to be released around the same time as PostGIS 3.0.

I think they are more or less sold on the idea except they did ask about
current market share and usage.

Part of the reason for that is our previous editions focused a lot on  "How
do I use this function or do this weird sounding thing that only GIS people
can make sense of"  instead of "How do I do this real world thing"

So one of the thoughts was having our table of contents be more like "How do
I do this with PostGIS" in somewhat laymen terms that most people can relate
to - like Political Districting, Real Estate analysis (walk scores,
elevation measurements to determine viablility of building on a plot of
land) 
 without scaring people off with "real world things" they can't relate to or
in overly techy terms.

Also since the 2nd Edition (which was in 2015 super ancient now since the
New shiny version at the time was 2.1 and 2.1 is not even supported
anymore).
Other major thing changed is a lot of people are deploying PostGIS on cloud
offerings like Amazon RDS, Microsoft Azure, and Google PostgreSQL for Cloud
so we plan to cover a bit about some things relevant in those that may not
be relevant when deploying on your own server.
 
That said, if people can respond with what things they are currently using
PostGIS for and also what hosting they are using for PostGIS, that would be
helpful for us to get a better idea of focus points.

It'd be great if you posted on the list, but if you are shy or need your
usage anonymized, you can write directly to me.

Thanks,
Regina

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

Re: [postgis-users] Split polygon by line

[Birgit Laggner]

Hi, 

I believe QGIS can not deal with GeometryCollections which are the resulting 
geometrytype of ST_Split. 

Your idea is not bad, but why do you save the result as wkt? 

I would try: 
SELECT (ST_Dump(ST_Split(circle, line))).geom::geometry (polygon,0) As geom 
INTO layerq 
FROM (SELECT 
ST_MakeLine(ST_MakePoint(10, 10),ST_MakePoint(190, 190)) As line, 
ST_Buffer(ST_GeomFromText('POINT(100 90)'), 50) As circle) As foo; 

If you cast the geometry type like that (::geometry (polygon,0)), PostGIS 
should be able to populate the geometry_columns table automaticly. 

Regards, 
Birgit 


Von: "Shane Carey"  
An: "arnaud listes"  
CC: postgis-users@lists.osgeo.org 
Gesendet: Donnerstag, 18. Oktober 2018 18:55:44 
Betreff: Re: [postgis-users] Split polygon by line 

Ok thanks, will do in future: 

This is what I have tried and no joy- any ideas? Thanks in advance. 

SELECT ST_AsText((ST_Dump(ST_Split(circle, line))).geom) As wkt 
INTO layerq 
FROM (SELECT 
ST_MakeLine(ST_MakePoint(10, 10),ST_MakePoint(190, 190)) As line, 
ST_Buffer(ST_GeomFromText('POINT(100 90)'), 50) As circle) As foo; 
select populate_geometry_columns(); 

Le gach dea ghui, 
Shane Carey 
GIS and Data Solutions Consultant 


On Thu, Oct 18, 2018 at 3:38 PM Arnaud L. < [ mailto:arnaud.lis...@codata.eu | 
arnaud.lis...@codata.eu ] > wrote: 


Le 18/10/2018 à 18:43, Shane Carey a écrit : 
> Yep, I run the following: 
> 
> SELECT ST_Split(circle, line) 
> INTO qlayer 
> FROM (SELECT 
> ST_MakeLine(ST_MakePoint(10, 10),ST_MakePoint(190, 190)) As line, 
> ST_Buffer(ST_GeomFromText('POINT(100 90)'), 50) As circle) As foo; 
> 
> But this layer does not show up as a gis layer in qgis? 

Probably because the geometry_columns table is not populated. 
Try to either create the destination table first and then to insert in 
it, or maybe just run : 
SELECT Populate_Geometry_Columns(); 

> Thanks in advance for your help. Sorry, what do you mean by top post? 

[ https://en.wikipedia.org/wiki/Posting_style | 
https://en.wikipedia.org/wiki/Posting_style ] 
Bottom or interleaved posting are the preferred style. 

-- 
Arnaud 




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

Re: [postgis-users] GEOSUnaryUnion: TopologyException: found non-noded intersection BUT geom is valid

[Birgit Laggner]

Dear Darafei, 

in preparation of filing the ticket, I tested the ST_Union on my selected, 
supposedly problematic geometries. Interestingly, the ST_Union worked without a 
problem on my small selection table. I aborted the ticket creation because I 
can't file this as a reproducible bug. 

My (then unnecessarily complicated workaround) of course worked fine, too. That 
means, I am able to get my work done, unfortunately, I still don't know what 
the original problem was. 

I am willing to let it rest and mark it as mysterious behaviour of a large 
dataset, as long as nobody has a specific interest to go to the bottom of the 
problem? 

Thanks a lot again for the idea with the wrapper to extract to problematic 
rows. 

Kind regards, 
Birgit 


Von: "Birgit Laggner"  
An: "PostGIS Users Discussion"  
Gesendet: Mittwoch, 19. September 2018 09:47:41 
Betreff: Re: [postgis-users] GEOSUnaryUnion: TopologyException: found non-noded 
intersection BUT geom is valid 

Dear Darafei, 

good idea to write a wrapper to ferret out the offending rows. I had to change 
it a bit because I don't use ST_Union(geometry, geometry) but the aggregate 
ST_Union(geometry). Since I am not very experienced in writing aggregates, I 
modified the wrapper that it fits for ST_Union(geometry[]) and used it in my 
query with ST_Union_Fails(array_agg(geom)). 

I have found 3 groups with a total of 923 geometries which returned an error 
with ST_Union. 

As soon as I get access to my old userID, I will file the them as a ticket. 

In the meantime I will try to union the problematic geometries in a loop one at 
a time with st_makevalid after each loop - maybe that will work as a 
workaround, now I could reduce the number of geometries dramatically. 


Thank you very much for your help! 

Kind regards, 
Birgit 




Von: "Darafei \"Komяpa\" Praliaskouski"  
An: "PostGIS Users Discussion"  
Gesendet: Mittwoch, 19. September 2018 02:46:42 
Betreff: Re: [postgis-users] GEOSUnaryUnion: TopologyException: found non-noded 
intersection BUT geom is valid 

Hi, 
This is not supposed to happen. 

Can you please isolate the offending rows unioning them pairwise and file them 
as a ticket? Here's a handy helper function that returns True if Union failed: 

create or replace function st_union_fails (geom1 geometry , geom2 geometry ) 
returns boolean 
as $$ 
begin 
begin 
geom1 = ST_Union (Geom1 , geom2) ; 
return false ; 
exception when others 
then 
return true ; 
end ; 
end ; 
$$ 
language plpgsql ; 



вт, 18 сент. 2018 г. в 16:25, Birgit Laggner < [ 
mailto:birgit.lagg...@thuenen.de | birgit.lagg...@thuenen.de ] >: 



Dear discussion group, 

I am trying to ST_Union several polygons and get a TopologyException 
(GEOSUnaryUnion: TopologyException: found non-noded intersection between 
LINESTRING (3.56442e+06 5.42679e+06, 3.56442e+06 5.42679e+06) and LINESTRING 
(3.56442e+06 5.42679e+06, 3.56442e+06 5.42679e+06) at 3564420.7944701263 
5426786.9800475985). 
I have tested with ST_Valid and all polygons in the datasets are valid. 

The dataset is pretty large (53.5 Mio. polygons) and will probably get grouped 
into around 8.7 Mio resulting geometries during ST_Union. 

If I try to ST_Union the geometries near the point ST_Union mentions in the 
error message (ST_Buffer with up to 100 metres), everything works fine. I am 
out of ideas how to find the problematic geometry or at least how to work 
around the TopologyException. 

My PostGIS version is: "POSTGIS="2.3.3 r15473" GEOS="3.5.1-CAPI-1.9.1 r4246" 
SFCGAL="1.2.2" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 
2015/09/16" LIBXML="2.9.3" LIBJSON="0.11.99" TOPOLOGY RASTER" 


I would be very glad if anyone would come up with ideas how to solve my 
problem. 

Thanks a lot in advance! 

Regards, 
Birgit 


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



-- 
Darafei Praliaskouski 
Support me: [ http://patreon.com/komzpa | http://patreon.com/komzpa ] 

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

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

Re: [postgis-users] GEOSUnaryUnion: TopologyException: found non-noded intersection BUT geom is valid

[Birgit Laggner]

Dear Darafei, 

good idea to write a wrapper to ferret out the offending rows. I had to change 
it a bit because I don't use ST_Union(geometry, geometry) but the aggregate 
ST_Union(geometry). Since I am not very experienced in writing aggregates, I 
modified the wrapper that it fits for ST_Union(geometry[]) and used it in my 
query with ST_Union_Fails(array_agg(geom)). 

I have found 3 groups with a total of 923 geometries which returned an error 
with ST_Union. 

As soon as I get access to my old userID, I will file the them as a ticket. 

In the meantime I will try to union the problematic geometries in a loop one at 
a time with st_makevalid after each loop - maybe that will work as a 
workaround, now I could reduce the number of geometries dramatically. 


Thank you very much for your help! 

Kind regards, 
Birgit 




Von: "Darafei \"Komяpa\" Praliaskouski"  
An: "PostGIS Users Discussion"  
Gesendet: Mittwoch, 19. September 2018 02:46:42 
Betreff: Re: [postgis-users] GEOSUnaryUnion: TopologyException: found non-noded 
intersection BUT geom is valid 

Hi, 
This is not supposed to happen. 

Can you please isolate the offending rows unioning them pairwise and file them 
as a ticket? Here's a handy helper function that returns True if Union failed: 

create or replace function st_union_fails (geom1 geometry , geom2 geometry ) 
returns boolean 
as $$ 
begin 
begin 
geom1 = ST_Union (Geom1 , geom2) ; 
return false ; 
exception when others 
then 
return true ; 
end ; 
end ; 
$$ 
language plpgsql ; 



вт, 18 сент. 2018 г. в 16:25, Birgit Laggner < [ 
mailto:birgit.lagg...@thuenen.de | birgit.lagg...@thuenen.de ] >: 



Dear discussion group, 

I am trying to ST_Union several polygons and get a TopologyException 
(GEOSUnaryUnion: TopologyException: found non-noded intersection between 
LINESTRING (3.56442e+06 5.42679e+06, 3.56442e+06 5.42679e+06) and LINESTRING 
(3.56442e+06 5.42679e+06, 3.56442e+06 5.42679e+06) at 3564420.7944701263 
5426786.9800475985). 
I have tested with ST_Valid and all polygons in the datasets are valid. 

The dataset is pretty large (53.5 Mio. polygons) and will probably get grouped 
into around 8.7 Mio resulting geometries during ST_Union. 

If I try to ST_Union the geometries near the point ST_Union mentions in the 
error message (ST_Buffer with up to 100 metres), everything works fine. I am 
out of ideas how to find the problematic geometry or at least how to work 
around the TopologyException. 

My PostGIS version is: "POSTGIS="2.3.3 r15473" GEOS="3.5.1-CAPI-1.9.1 r4246" 
SFCGAL="1.2.2" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 
2015/09/16" LIBXML="2.9.3" LIBJSON="0.11.99" TOPOLOGY RASTER" 


I would be very glad if anyone would come up with ideas how to solve my 
problem. 

Thanks a lot in advance! 

Regards, 
Birgit 


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



-- 
Darafei Praliaskouski 
Support me: [ http://patreon.com/komzpa | http://patreon.com/komzpa ] 

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

[postgis-users] GEOSUnaryUnion: TopologyException: found non-noded intersection BUT geom is valid

[Birgit Laggner]

Dear discussion group, 

I am trying to ST_Union several polygons and get a TopologyException 
(GEOSUnaryUnion: TopologyException: found non-noded intersection between 
LINESTRING (3.56442e+06 5.42679e+06, 3.56442e+06 5.42679e+06) and LINESTRING 
(3.56442e+06 5.42679e+06, 3.56442e+06 5.42679e+06) at 3564420.7944701263 
5426786.9800475985). 
I have tested with ST_Valid and all polygons in the datasets are valid. 

The dataset is pretty large (53.5 Mio. polygons) and will probably get grouped 
into around 8.7 Mio resulting geometries during ST_Union. 

If I try to ST_Union the geometries near the point ST_Union mentions in the 
error message (ST_Buffer with up to 100 metres), everything works fine. I am 
out of ideas how to find the problematic geometry or at least how to work 
around the TopologyException. 

My PostGIS version is: "POSTGIS="2.3.3 r15473" GEOS="3.5.1-CAPI-1.9.1 r4246" 
SFCGAL="1.2.2" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 
2015/09/16" LIBXML="2.9.3" LIBJSON="0.11.99" TOPOLOGY RASTER" 


I would be very glad if anyone would come up with ideas how to solve my 
problem. 

Thanks a lot in advance! 

Regards, 
Birgit 


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

Re: [postgis-users] MultiPolygon -> Polygon

[Birgit Laggner]

Hi James, 

perhaps you could split the data per multipolygon into the largest polygon and 
the rest. Then buffer the smaller polygons by a small amount and union them 
with the large polygon? 

Regards, 
Birgit 


Von: "James Sewell"  
An: postgis-users@lists.osgeo.org 
Gesendet: Montag, 6. August 2018 09:24:42 
Betreff: [postgis-users] MultiPolygon -> Polygon 

Hi all, 
I have a set of LineStrings which are being turned in Polygons (most of which 
are invalid). 

I can clean these (makevalid or buffer), but then some of them result in 
MultiPolygons (which makes sense - mostly small bowties). 

I'm looking for the easiest way to get a single Polygon from each of these 
MultiPolygons - I can tolerate a small amount of change in area (up to 5%). 

I've tried small just keeping the largest poly, buffer(small_number) and 
concavehull but nothing seems quite right. 

Any ideas? 

Cheers, 

James Sewell, 




The contents of this email are confidential and may be subject to legal or 
professional privilege and copyright. No representation is made that this email 
is free of viruses or other defects. If you have received this communication in 
error, you may not copy or distribute any part of it or otherwise disclose its 
contents to anyone. Please advise the sender of your incorrect receipt of this 
correspondence. 




The contents of this email are confidential and may be subject to legal or 
professional privilege and copyright. No representation is made that this email 
is free of viruses or other defects. If you have received this communication in 
error, you may not copy or distribute any part of it or otherwise disclose its 
contents to anyone. Please advise the sender of your incorrect receipt of this 
correspondence. 



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

Re: [postgis-users] ST_Union results in empty polygon

[Birgit Laggner]

Hi Jonathan, 

I tested your example and it seems the reason why the ST_Union results in an 
empty polygon is because the resulting polygons of the ST_Buffer are already 
empty. I could imagine complex polygons don't go well with a very small buffer. 

In my experience, ST_Union only results in empty polygons when the input 
polygons are not valid. Did you test if your polygons are valid? If validity is 
indeed the problem, you can solve the issue with ST_MakeValid rather than with 
ST_Buffer. 

Regards, 

Birgit 


Von: "Jonathan McCormack"  
An: postgis-users@lists.osgeo.org 
Gesendet: Montag, 26. März 2018 20:05:42 
Betreff: [postgis-users] ST_Union results in empty polygon 



We've seen occasional inconsistencies from ST_Union where the resulting polygon 
is empty or clearly cannot be the result of the constituent geometries. These 
issues can often be resolved by applying a very slight buffer, though sometimes 
the reverse is true, where the ST_Union of the unbuffered geometries is fine, 
but applying the buffer causes the incorrect geometry. The following query 
shows this latter case (pardon the lengthy geometry strings): 



WITH poly AS ( 

SELECT 

ST_SetSRID( 

Re: [postgis-users] raster query

[Birgit Laggner]

I am glad I could help :-)

Regards, Birgit

Am 04.11.2016 um 15:45 schrieb Stephen Crawford:
I had to make just a few minor adjustments, only about casting types 
which you couldn't have know about anyway.  Worked well. Many Thanks.


On 11/3/2016 8:20 AM, Birgit Laggner wrote:

Hi Stephen,

here my proposed solution (untested!!)

DO $$
DECLARE r record; doy date;
BEGIN

FOR r IN (SELECT DISTINCT ST_Value(rast,x,y) AS obs_date FROM 
doy_raster, generate_series(1,(SELECT ST_Width(rast) FROM 
doy_raster)) AS x, generate_series(1,(SELECT ST_Height(rast) FROM 
doy_raster)) AS y) LOOP


doy := r.obs_date;

EXECUTE
'INSERT INTO accum_risk (obs_date, rast)
 SELECT '||doy||',
ST_Union(ST_MapAlgebra(a.rast, b.rast,
 ''CASE WHEN [rast2]='||quote_literal(doy)||' 
THEN [rast1] ELSE 0 END'',

 ''64BF'',''INTERSECTION''),''SUM'')
 FROM daily_risk AS a, doy_raster AS b
 WHERE a.obs_date BETWEEN '||quote_literal(doy)||'::date - ''60 
days''::interval AND

  '||quote_literal(doy);

END LOOP;

END $$;

The idea is to loop over the all existing observation dates from 
doy_raster and use ST_MapAlgebra to limit the raster cells from 
daily_risk raster to those intersecting with raster cells from 
doy_raster with the corresponding observation date for the current loop.
As already mentioned, the code is not tested (because I don't have 
similar datasets available)...


Regards,

Birgit


Am 02.11.2016 um 21:51 schrieb Stephen Crawford:

Hi All,

I am hoping somebody can help me with a query.  I have a table of 
rasters where each record is date ("obs_date") and raster containing 
a risk value of 0  or 1.  My easy, successful query to accumulate 
the risk values over the previous 60 days is:


INSERT INTO accum_risk (obs_date, rast)
SELECT '1979-07-15', ST_Union(rast,'SUM')
FROM daily_risk
WHERE obs_date BETWEEN '1979-05-16' AND '1979-07-15';

My final goal--for which I am asking help--is similar to the above 
query, but it will reference another raster table.  This table has 
for each grid cell a value for the day of year (DOY) from which the 
60 day accumulation should be made. Conceptually:


FOR EACH doy_cell IN doy_raster
SELECT ST_Union(rast,'SUM')
FROM daily_risk
WHERE obs_date BETWEEN doy-60 AND doy;

Any help is greatly appreciated.

Thanks,
Steve
--
Stephen Crawford
Center for Environmental Informatics
The Pennsylvania State University
src...@psu.edu


___
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


--
Stephen Crawford
Center for Environmental Informatics
The Pennsylvania State University
src...@psu.edu
814.865.9905


___
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] raster query

[Birgit Laggner]

Hi Stephen,

here my proposed solution (untested!!)

DO $$
DECLARE r record; doy date;
BEGIN

FOR r IN (SELECT DISTINCT ST_Value(rast,x,y) AS obs_date FROM 
doy_raster, generate_series(1,(SELECT ST_Width(rast) FROM doy_raster)) 
AS x, generate_series(1,(SELECT ST_Height(rast) FROM doy_raster)) AS y) LOOP


doy := r.obs_date;

EXECUTE
'INSERT INTO accum_risk (obs_date, rast)
 SELECT '||doy||',
ST_Union(ST_MapAlgebra(a.rast, b.rast,
 ''CASE WHEN [rast2]='||quote_literal(doy)||' THEN 
[rast1] ELSE 0 END'',

 ''64BF'',''INTERSECTION''),''SUM'')
 FROM daily_risk AS a, doy_raster AS b
 WHERE a.obs_date BETWEEN '||quote_literal(doy)||'::date - ''60 
days''::interval AND

  '||quote_literal(doy);

END LOOP;

END $$;

The idea is to loop over the all existing observation dates from 
doy_raster and use ST_MapAlgebra to limit the raster cells from 
daily_risk raster to those intersecting with raster cells from 
doy_raster with the corresponding observation date for the current loop.
As already mentioned, the code is not tested (because I don't have 
similar datasets available)...


Regards,

Birgit


Am 02.11.2016 um 21:51 schrieb Stephen Crawford:

Hi All,

I am hoping somebody can help me with a query.  I have a table of 
rasters where each record is date ("obs_date") and raster containing a 
risk value of 0  or 1.  My easy, successful query to accumulate the 
risk values over the previous 60 days is:


INSERT INTO accum_risk (obs_date, rast)
SELECT '1979-07-15', ST_Union(rast,'SUM')
FROM daily_risk
WHERE obs_date BETWEEN '1979-05-16' AND '1979-07-15';

My final goal--for which I am asking help--is similar to the above 
query, but it will reference another raster table.  This table has for 
each grid cell a value for the day of year (DOY) from which the 60 day 
accumulation should be made. Conceptually:


FOR EACH doy_cell IN doy_raster
SELECT ST_Union(rast,'SUM')
FROM daily_risk
WHERE obs_date BETWEEN doy-60 AND doy;

Any help is greatly appreciated.

Thanks,
Steve
--
Stephen Crawford
Center for Environmental Informatics
The Pennsylvania State University
src...@psu.edu


___
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

[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 
<birgit.lagg...@thuenen.de <mailto:birgit.lagg...@thuenen.de>> 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
<godzilalal...@gmail.com <mailto:godzilalal...@gmail.com>> 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 <http://line.id>;

END LOOP;

END

$$;


>
> On Tue, Sep 13, 2016 at 11:12 AM James Keener
<j...@jimkeener.com <mailto:j...@jimkeener.com>> 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
<mailto:jon.mala...@shrewsburyma.gov>> wrote:
>
> 

Re: [postgis-users] st_clusterwithin on postgis 2

[Birgit Laggner]

Hi Paolo,

sorry, I did not really think it through, but now the changed and tested 
example:


DO $$
DECLARE distance numeric; t_schema varchar; input_table varchar; 
output_table varchar; r record; gid integer; geom geometry; c record; 
c_gid integer[]; c_geom geometry[]; i record;

BEGIN

distance := 200;
t_schema := 'p_sam';
input_table := 'cluster_test_input';
output_table := 'cluster_test_output';


CREATE TEMP TABLE clustered (gid integer) ON COMMIT DROP;

EXECUTE 'SELECT gid, geom FROM 
'||quote_ident(t_schema)||'.'||quote_ident(input_table) INTO r;


gid := r.gid;
geom := r.geom;

EXECUTE 'INSERT INTO clustered SELECT '||gid;

LOOP

EXECUTE 'SELECT array_agg(gid) as gid, array_agg(geom) as geom
 FROM 
'||quote_ident(t_schema)||'.'||quote_ident(input_table)||'
 WHERE gid NOT IN (SELECT gid FROM clustered) AND 
ST_DWithin($1, geom, $2)'

 INTO c USING geom, distance;

c_gid := c.gid;
c_geom := c.geom || array[geom]::geometry[];

EXECUTE 'INSERT INTO clustered SELECT unnest($1)' USING c_gid;

geom := ST_ForceCollection(ST_Collect(c_geom));

EXIT WHEN c IS NULL;

END LOOP;

EXECUTE 'INSERT INTO 
'||quote_ident(t_schema)||'.'||quote_ident(output_table)||' VALUES($1)' 
USING geom;


LOOP
BEGIN
EXECUTE 'SELECT gid, geom FROM 
'||quote_ident(t_schema)||'.'||quote_ident(input_table)||'

 WHERE gid NOT IN (SELECT gid FROM clustered)' INTO r;

gid := r.gid;
geom := r.geom;

EXECUTE 'INSERT INTO clustered SELECT '||gid;

LOOP

EXECUTE 'SELECT array_agg(gid) as gid, array_agg(geom) as geom
 FROM 
'||quote_ident(t_schema)||'.'||quote_ident(input_table)||'
 WHERE gid NOT IN (SELECT gid FROM clustered) AND 
ST_DWithin($1, geom, $2)'

 INTO c USING geom, distance;

c_gid := c.gid;
c_geom := c.geom || array[geom]::geometry[];

EXECUTE 'INSERT INTO clustered SELECT unnest($1)' USING c_gid;

geom := ST_ForceCollection(ST_Collect(c_geom));

EXIT WHEN c IS NULL;

END LOOP;

EXECUTE 'INSERT INTO 
'||quote_ident(t_schema)||'.'||quote_ident(output_table)||'

 VALUES ($1)' USING geom;

EXCEPTION WHEN SQLSTATE '22004' THEN EXIT;
END;
END LOOP;

END $$


Regards,

Birgit



Am 21.09.2016 um 14:53 schrieb Birgit Laggner:

Hi Paolo,

here an example how it could work (mind that you have to replace the 
??). You could also recode the DO block into a function if you would 
like that better.


DO $$
DECLARE distance numeric; t_schema varchar; input_table varchar; 
output_table varchar; r record; gid integer; geom geometry; c record; 
c_gid integer[]; c_geom geometry;

BEGIN

distance := ??;
t_schema := '??';
input_table := '??';
output_table := '??';


CREATE TEMP TABLE clustered (gid integer) ON COMMIT DROP;

EXECUTE 'SELECT gid, geom FROM 
'||quote_ident(t_schema)||'.'||quote_ident(input_table) INTO r;


gid := r.gid;
geom := r.geom;

EXECUTE 'INSERT INTO clustered SELECT '||gid;

EXECUTE 'SELECT array_agg(gid) as gid, 
ST_ForceCollection(ST_Collect(ST_Collect(geom),$2)) as geom FROM 
'||quote_ident(t_schema)||'.'||quote_ident(input_table)||' WHERE gid 
!= $1 AND ST_DWithin($2, geom, $3)' INTO c USING gid, geom, distance;


c_gid := c.gid;
c_geom := c.geom;

EXECUTE 'INSERT INTO 
'||quote_ident(t_schema)||'.'||quote_ident(output_table)||' 
VALUES($1)' USING c_geom;


EXECUTE 'INSERT INTO clustered SELECT unnest($1)' USING c_gid;

LOOP
BEGIN
EXECUTE 'SELECT gid, geom FROM 
'||quote_ident(t_schema)||'.'||quote_ident(input_table)||'

 WHERE gid NOT IN (SELECT gid FROM clustered)' INTO r;

gid := r.gid;
geom := r.geom;

EXECUTE 'INSERT INTO clustered SELECT '||gid;

EXECUTE 'SELECT array_agg(gid) as gid, 
ST_ForceCollection(ST_Collect(ST_Collect(geom),$1)) as geom

 FROM '||quote_ident(t_schema)||'.'||quote_ident(input_table)||'
 WHERE gid NOT IN (SELECT gid FROM clustered) AND 
ST_DWithin($1, geom, $2)'

 INTO c USING geom, distance;

c_gid := c.gid;
c_geom := c.geom;

EXECUTE 'INSERT INTO 
'||quote_ident(t_schema)||'.'||quote_ident(output_table)||'

 VALUES ($1)' USING c_geom;

EXECUTE 'INSERT INTO clustered SELECT unnest($1)' USING c_gid;

EXCEPTION WHEN SQLSTATE '22004' THEN EXIT;
END;
END LOOP;

END $$

Regards,

Birgit
Am 21.09.2016 um 07:47 schrieb Paolo Importuni:

Hi all,
I need to run a query on my postgres/postgis (9.3/2.0) that uses 
ST_CLUSTERWITHIN function. This function is available since postgis 
2.2.0 so I should do a soft upgrade on my ubuntu server.  Since we 
have a kind of customer demo in  a few days and have no much time 
left, I am not willing to change our setup right now, and I'd rather 
do it after that meeting.
The question is: is there a way to aggregate geometries like 
ST_ClusterWIthin does? Anybody can provide any working examples?


Thanks and regards

Paolo I

Re: [postgis-users] st_clusterwithin on postgis 2

[Birgit Laggner]

Hi Paolo,

here an example how it could work (mind that you have to replace the 
??). You could also recode the DO block into a function if you would 
like that better.


DO $$
DECLARE distance numeric; t_schema varchar; input_table varchar; 
output_table varchar; r record; gid integer; geom geometry; c record; 
c_gid integer[]; c_geom geometry;

BEGIN

distance := ??;
t_schema := '??';
input_table := '??';
output_table := '??';


CREATE TEMP TABLE clustered (gid integer) ON COMMIT DROP;

EXECUTE 'SELECT gid, geom FROM 
'||quote_ident(t_schema)||'.'||quote_ident(input_table) INTO r;


gid := r.gid;
geom := r.geom;

EXECUTE 'INSERT INTO clustered SELECT '||gid;

EXECUTE 'SELECT array_agg(gid) as gid, 
ST_ForceCollection(ST_Collect(ST_Collect(geom),$2)) as geom FROM 
'||quote_ident(t_schema)||'.'||quote_ident(input_table)||' WHERE gid != 
$1 AND ST_DWithin($2, geom, $3)' INTO c USING gid, geom, distance;


c_gid := c.gid;
c_geom := c.geom;

EXECUTE 'INSERT INTO 
'||quote_ident(t_schema)||'.'||quote_ident(output_table)||' VALUES($1)' 
USING c_geom;


EXECUTE 'INSERT INTO clustered SELECT unnest($1)' USING c_gid;

LOOP
BEGIN
EXECUTE 'SELECT gid, geom FROM 
'||quote_ident(t_schema)||'.'||quote_ident(input_table)||'

 WHERE gid NOT IN (SELECT gid FROM clustered)' INTO r;

gid := r.gid;
geom := r.geom;

EXECUTE 'INSERT INTO clustered SELECT '||gid;

EXECUTE 'SELECT array_agg(gid) as gid, 
ST_ForceCollection(ST_Collect(ST_Collect(geom),$1)) as geom

 FROM '||quote_ident(t_schema)||'.'||quote_ident(input_table)||'
 WHERE gid NOT IN (SELECT gid FROM clustered) AND 
ST_DWithin($1, geom, $2)'

 INTO c USING geom, distance;

c_gid := c.gid;
c_geom := c.geom;

EXECUTE 'INSERT INTO 
'||quote_ident(t_schema)||'.'||quote_ident(output_table)||'

 VALUES ($1)' USING c_geom;

EXECUTE 'INSERT INTO clustered SELECT unnest($1)' USING c_gid;

EXCEPTION WHEN SQLSTATE '22004' THEN EXIT;
END;
END LOOP;

END $$

Regards,

Birgit

Am 21.09.2016 um 07:47 schrieb Paolo Importuni:

Hi all,
I need to run a query on my postgres/postgis (9.3/2.0) that uses 
ST_CLUSTERWITHIN function. This function is available since postgis 
2.2.0 so I should do a soft upgrade on my ubuntu server.  Since we 
have a kind of customer demo in  a few days and have no much time 
left, I am not willing to change our setup right now, and I'd rather 
do it after that meeting.
The question is: is there a way to aggregate geometries like 
ST_ClusterWIthin does? Anybody can provide any working examples?


Thanks and regards

Paolo I.


___
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

[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] ST_Split with Multilinestring

[Birgit Laggner]

Hi Sandro,

mmh... I really would like to, but I am no programmer and wouldn't know 
how to contribute and I am not in a position to decide about financing 
somebody to add this functionality. I'm sorry.


But thank you anyway for clarifying the current functionality status!

Regards,

Birgit

Am 16.06.2016 um 10:34 schrieb Sandro Santilli:

On Wed, Jun 15, 2016 at 07:19:34PM +0200, Birgit Laggner, vTI wrote:

Hi Sandro,

yes, I read so, too. But, our database has PostGIS version:
"POSTGIS="2.2.1 r14555" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel.
4.9.1, 04 March 2015" GDAL="GDAL 1.11.2, released 2015/02/10"
LIBXML="2.9.2" LIBJSON="0.11.99" TOPOLOGY RASTER" and I still get
that error. Any idea why?

Oops, sorry: 2.2.0 added capability to split _lines_ by multilines,
not _polygons_. Polygons are still only splittable by single linestrings.

It should be an easy addition, if you want to contribute or finance one.

--strk;



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

Re: [postgis-users] ST_Split with Multilinestring

[Birgit Laggner]

Hi Rémi,

thank you for sharing your function. In the meantime I seem to have 
solved the problem by combining the st_split function with a function 
from an old use case:


CREATE OR REPLACE FUNCTION ST_SplitPolygon(poly geometry, blade 
geometry) RETURNS geometry AS

$$
declare
out_geom geometry;
begin
begin
out_geom := st_split(poly, blade);
exception when internal_error then
out_geom := ST_Collect((d).geom) FROM (SELECT 
ST_Dump(ST_Polygonize(ST_Union(ST_Boundary(poly), blade))) AS d) f WHERE 
ST_Area(ST_Intersection(poly, (d).geom)) / ST_Area((d).geom) >= 1 - 1e-10;

end;
return out_geom;
end;
$$
LANGUAGE plpgsql IMMUTABLE STRICT COST 100;

But I will definitely give your function a try and compare the results!

Thanks again! Regards,

Birgit

Am 16.06.2016 um 10:18 schrieb Rémi Cura:

Hey,
I wrote a multi-friendly st_split a while ago.
Maybe it will help you:
https://github.com/Remi-C/_utilities/blob/master/postgis/rc_split_multi.sql

Cheers,
Rémi-C

2016-06-15 19:19 GMT+02:00 Birgit Laggner, vTI 
<birgit.lagg...@thuenen.de <mailto:birgit.lagg...@thuenen.de>>:


Hi Sandro,

yes, I read so, too. But, our database has PostGIS version:
"POSTGIS="2.2.1 r14555" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel.
4.9.1, 04 March 2015" GDAL="GDAL 1.11.2, released 2015/02/10"
LIBXML="2.9.2" LIBJSON="0.11.99" TOPOLOGY RASTER" and I still get
that error. Any idea why?

Regards,

Birgit


Am 15.06.2016 um 19:02 schrieb Sandro Santilli:

On Tue, Jun 14, 2016 at 02:52:52PM +0200, Birgit Laggner wrote:

ERROR: Splitting a Polygon by a MultiLineString is unsupported

Support for this functionality was added in PostGIS-2.2.0

--strk;



___
postgis-users mailing list
postgis-users@lists.osgeo.org <mailto: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] ST_Split with Multilinestring

[Birgit Laggner, vTI]

Hi Sandro,

yes, I read so, too. But, our database has PostGIS version: 
"POSTGIS="2.2.1 r14555" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.9.1, 
04 March 2015" GDAL="GDAL 1.11.2, released 2015/02/10" LIBXML="2.9.2" 
LIBJSON="0.11.99" TOPOLOGY RASTER" and I still get that error. Any idea why?


Regards,

Birgit

Am 15.06.2016 um 19:02 schrieb Sandro Santilli:

On Tue, Jun 14, 2016 at 02:52:52PM +0200, Birgit Laggner wrote:


ERROR: Splitting a Polygon by a MultiLineString is unsupported

Support for this functionality was added in PostGIS-2.2.0

--strk;




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

Re: [postgis-users] ST_Split with Multilinestring

[Birgit Laggner]

Yes, this function will probably work - that's similar to what I meant 
with "looping the ST_Split with some function, DO block or pgScript 
after selecting all intersecting linestrings per polygon" in my original 
mail. If my current approach fails, too, I will revert to some sort of 
looping approach like the one you suggested. Thanks again for taking the 
time to analyze my problem!


Regards,

Birgit


Achtung: Das Thünen-Institut hat die Domain gewechselt.
 Bitte ändern Sie meine Mailadresse in Ihrem Adressbuch!

Notice: Thünen Institute has changed its domain.
Please change my email address in your address book!


Dipl.-Geoökol. Birgit Laggner

Thünen-Institut für Ländliche Räume / Thünen Institute of Rural Studies
Bundesallee 50
38116 Braunschweig (Germany)

Tel: +49 531 596-5240
Mail: birgit.lagg...@thuenen.de
Web: www.thuenen.de

Am 15.06.2016 um 11:27 schrieb Marcin Mionskowski:
I think I finaly understand what you have and what you want to achieve 
- take a look at polygon and multiline definitions :)


If I'm right, this function should do the work:

create table polygon
(
id serial,
geom geometry,
done smallint default 0
);
create table multiline
(
id serial,
geom geometry,
done smallint default 0
);
create table line
(
id serial,
geom geometry,
done smallint default 0
);

insert into polygon (geom)
select ST_Buffer(ST_GeomFromText('POINT(10 10)'), 10);

insert into multiline (geom)
select ST_Union(
(select ST_MakeLine(ST_MakePoint(0, 1),ST_MakePoint(20, 21))),
(select ST_MakeLine(ST_MakePoint(1, 0),ST_MakePoint(20, 19)))
);

create or replace function test() returns void as
$BODY$

declare ln geometry;
i int;

BEGIN
truncate table line;

insert into line (geom)
select (ST_Dump(geom)).geom
from multiline;

loop
select geom,id
into ln,i
from line
where done=0
limit 1;

with
del as (
delete from polygon p
where ST_Intersects(p.geom,ln)
returning p.geom
)
insert into polygon (geom)
select (st_dump(st_split(del.geom,ln))).geom
from del;

update line
set done=1
where id=i;

exit when i is null;
end loop;
end;
$BODY$
LANGUAGE plpgsql;

select test();



On Wed, 15 Jun 2016 11:10:33 +0200, Birgit Laggner 
<birgit.lagg...@thuenen.de> wrote:



Dear list,

@Marcone: Unfortunately, the st_multi didn't change anything. The error
message is still the same.

Next, I will try some function from a user example from 2008
(https://trac.osgeo.org/postgis/wiki/UsersWikiSplitPolygonWithLineString). 


Maybe this will be a workaround.

But I am still wondering why I get the error message at all. I did find
a PostGIS ticket addressing the same functionality
(https://trac.osgeo.org/postgis/ticket/3097) and there is mentioned that
the problem is fixed with PostGIS 2.2.0. My PostGIS version is
"POSTGIS="2.2.1 r14555" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.9.1,
04 March 2015" GDAL="GDAL 1.11.2, released 2015/02/10" LIBXML="2.9.2"
LIBJSON="0.11.99" TOPOLOGY RASTER". So, I would suppose I shouldn't have
this problem.

Any other ideas?

Thanks and regards,

Birgit

Am 15.06.2016 um 10:01 schrieb Birgit Laggner:

Hi Marcin,

thanks for your reply. I don't think, that will provide the results I
am looking for. I am deliberately aggregating all intersecting lines
with ST_Union because otherwise the query would result in several
collections per polygon, each collection containig the st_split result
for one of the intersecting lines (blades). But I want the polygon cut
with all blades.

Currently I am trying Marcone's suggestion (although I cannot imagine
why this should solve anything - but you never know until you try,
right?).

So, thank you both, Marcin and Marcone!

Regards,

Birgit

Am 15.06.2016 um 09:24 schrieb Marcin Mionskowski:

On Tue, 14 Jun 2016 15:37:40 +0200, Marcone <marconepe...@gmail.com>
wrote:



2016-06-14 9:52 GMT-03:00 Birgit Laggner <birgit.lagg...@thuenen.de>:

I would like to use the ST_Split function to split polygons from one
   table with all intersecting lines from another table.
Unfortunately,
   I get the following error:

 ERROR: Splitting a Polygon by a MultiLineString is unsupported

   SQL Status:XX000

 This is my query:

 select betr_id, (cut).path[1], (cut).geom from (select betr_id,
   st_dump(st_split) as cut from (select a.betr_id, ST_Split(a.geom,
   st_union(b.geom)) from p_sam.nihb_2013_convex_hull_betr a left 
join
   p_sam.ni_dlm13_aaa_gew_sie_ver_l b on st_intersects(a.geom, 
b.geom)

   group by a.betr_id, a.geom) sel1) sel2;

 The Manual contains the following info regarding this problem:

 "The function

Re: [postgis-users] ST_Split with Multilinestring

[Birgit Laggner]

Dear list,

@Marcone: Unfortunately, the st_multi didn't change anything. The error 
message is still the same.


Next, I will try some function from a user example from 2008 
(https://trac.osgeo.org/postgis/wiki/UsersWikiSplitPolygonWithLineString). 
Maybe this will be a workaround.


But I am still wondering why I get the error message at all. I did find 
a PostGIS ticket addressing the same functionality 
(https://trac.osgeo.org/postgis/ticket/3097) and there is mentioned that 
the problem is fixed with PostGIS 2.2.0. My PostGIS version is 
"POSTGIS="2.2.1 r14555" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.9.1, 
04 March 2015" GDAL="GDAL 1.11.2, released 2015/02/10" LIBXML="2.9.2" 
LIBJSON="0.11.99" TOPOLOGY RASTER". So, I would suppose I shouldn't have 
this problem.


Any other ideas?

Thanks and regards,

Birgit

Am 15.06.2016 um 10:01 schrieb Birgit Laggner:

Hi Marcin,

thanks for your reply. I don't think, that will provide the results I 
am looking for. I am deliberately aggregating all intersecting lines 
with ST_Union because otherwise the query would result in several 
collections per polygon, each collection containig the st_split result 
for one of the intersecting lines (blades). But I want the polygon cut 
with all blades.


Currently I am trying Marcone's suggestion (although I cannot imagine 
why this should solve anything - but you never know until you try, 
right?).


So, thank you both, Marcin and Marcone!

Regards,

Birgit

Am 15.06.2016 um 09:24 schrieb Marcin Mionskowski:
On Tue, 14 Jun 2016 15:37:40 +0200, Marcone <marconepe...@gmail.com> 
wrote:




2016-06-14 9:52 GMT-03:00 Birgit Laggner <birgit.lagg...@thuenen.de>:

I would like to use the ST_Split function to split polygons from one
   table with all intersecting lines from another table. 
Unfortunately,

   I get the following error:

 ERROR: Splitting a Polygon by a MultiLineString is unsupported

   SQL Status:XX000

 This is my query:

 select betr_id, (cut).path[1], (cut).geom from (select betr_id,
   st_dump(st_split) as cut from (select a.betr_id, ST_Split(a.geom,
   st_union(b.geom)) from p_sam.nihb_2013_convex_hull_betr a left join
   p_sam.ni_dlm13_aaa_gew_sie_ver_l b on st_intersects(a.geom, b.geom)
   group by a.betr_id, a.geom) sel1) sel2;

 The Manual contains the following info regarding this problem:

 "The function supports splitting a line by (multi)point, 
(multi)line

   or (multi)polygon boundary, a (multi)polygon by line."

 This might mean that I would be able to split a line by
   multilinestrings, but a polygon only by single linestrings - is 
that

   correct? Does anyone has a suggestion how I could work around this
   problem (aside from looping the ST_Split with some function, DO
   block or pgScript after selecting all intersecting linestrings per
   polygon)?

 Thanks a lot for any helpful suggestions!


I'm not sure if I understand your problem, but try use 
ST_Split(st_multi(a.geom),
   st_union(b.geom)). I'm not test this, but I think that will solve 
your problem.


Best regards.




Have you tried to ST_Dump multilinestrings first?
Using CTE this could look like:

with
a as ( --polygon
select ST_Buffer(ST_GeomFromText('POINT(100 90)'), 50) a
),
b as ( --multilinestring
select ST_Union(
(select ST_MakeLine(ST_MakePoint(10, 10),ST_MakePoint(190, 
190))),
(select ST_MakeLine(ST_MakePoint(191, 191),ST_MakePoint(200, 
200)))

) b
),
c as ( --multilinestring dump -> 2 separate linestrings
select (ST_Dump(b)).geom c
from b
)

SELECT ST_AsText(ST_Split(a, c))
from a,c
where ST_Intersects(a,c)

Intersection check is done here on line parts.

Regards
___
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] ST_Split with Multilinestring

[Birgit Laggner]

Hi Marcin,

thanks for your reply. I don't think, that will provide the results I am 
looking for. I am deliberately aggregating all intersecting lines with 
ST_Union because otherwise the query would result in several collections 
per polygon, each collection containig the st_split result for one of 
the intersecting lines (blades). But I want the polygon cut with all blades.


Currently I am trying Marcone's suggestion (although I cannot imagine 
why this should solve anything - but you never know until you try, right?).


So, thank you both, Marcin and Marcone!

Regards,

Birgit

Am 15.06.2016 um 09:24 schrieb Marcin Mionskowski:
On Tue, 14 Jun 2016 15:37:40 +0200, Marcone <marconepe...@gmail.com> 
wrote:




2016-06-14 9:52 GMT-03:00 Birgit Laggner <birgit.lagg...@thuenen.de>:

I would like to use the ST_Split function to split polygons from one
   table with all intersecting lines from another table. Unfortunately,
   I get the following error:

 ERROR: Splitting a Polygon by a MultiLineString is unsupported

   SQL Status:XX000

 This is my query:

 select betr_id, (cut).path[1], (cut).geom from (select betr_id,
   st_dump(st_split) as cut from (select a.betr_id, ST_Split(a.geom,
   st_union(b.geom)) from p_sam.nihb_2013_convex_hull_betr a left join
   p_sam.ni_dlm13_aaa_gew_sie_ver_l b on st_intersects(a.geom, b.geom)
   group by a.betr_id, a.geom) sel1) sel2;

 The Manual contains the following info regarding this problem:

 "The function supports splitting a line by (multi)point, 
(multi)line

   or (multi)polygon boundary, a (multi)polygon by line."

 This might mean that I would be able to split a line by
   multilinestrings, but a polygon only by single linestrings - is that
   correct? Does anyone has a suggestion how I could work around this
   problem (aside from looping the ST_Split with some function, DO
   block or pgScript after selecting all intersecting linestrings per
   polygon)?

 Thanks a lot for any helpful suggestions!


I'm not sure if I understand your problem, but try use 
ST_Split(st_multi(a.geom),
   st_union(b.geom)). I'm not test this, but I think that will solve 
your problem.


Best regards.




Have you tried to ST_Dump multilinestrings first?
Using CTE this could look like:

with
a as ( --polygon
select ST_Buffer(ST_GeomFromText('POINT(100 90)'), 50) a
),
b as ( --multilinestring
select ST_Union(
(select ST_MakeLine(ST_MakePoint(10, 10),ST_MakePoint(190, 
190))),
(select ST_MakeLine(ST_MakePoint(191, 191),ST_MakePoint(200, 
200)))

) b
),
c as ( --multilinestring dump -> 2 separate linestrings
select (ST_Dump(b)).geom c
from b
)

SELECT ST_AsText(ST_Split(a, c))
from a,c
where ST_Intersects(a,c)

Intersection check is done here on line parts.

Regards
___
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] Point In Polygon Update Query

[Birgit Laggner]

Hi Chris,

for the warning if a user is within a mile of the warned area, you could 
use ST_DWithin instead of ST_Within.


Since you have EPSG 4326, I would convert the geometries into geography 
types because then you can define the distance in meters (not in 
degrees), perhaps like:


ST_DWithin(mastergistmp.geom::geography, wwmaster.geom::geography, 1610.0)

Otherwise your UPDATE queries look fine for me...

Regards,

Birgit

Am 27.05.2016 um 02:04 schrieb Chris B:
I'm trying to put what the warning is for the point in the polygon, I 
have my wwmaster which is where all the warnings are stored, my 
mastergistmp where my users are stored, I know if a user is in a 
severe thunderstorm warning and a tornado warning the tornado warning 
will override the severe storm warning which is what I want does 
anyone see anything wrong with this I think I am missing some users in 
the warned area. everything is epsg 4326 I would like to also make it 
where users with in a mile of the warned area are updated with the 
information I've tried st_buffer and have had trouble with it. I'm 
updating the mastergistmp.warn or thats what I'm shooting for.


UPDATE mastergistmp SET warn = wwmaster.prod_type FROM wwmaster WHERE 
 wwmaster.prod_type = 'Flash Flood Warning' and 
ST_Within(mastergistmp.geom,wwmaster.geom);
UPDATE mastergistmp SET warn = wwmaster.prod_type FROM wwmaster WHERE 
 wwmaster.prod_type = 'Severe Thunderstorm Warning' and 
ST_Within(mastergistmp.geom,wwmaster.geom);
UPDATE mastergistmp SET warn = wwmaster.prod_type FROM wwmaster WHERE 
 wwmaster.prod_type = 'Tornado Warning' and 
ST_Within(mastergistmp.geom,wwmaster.geom);




___
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] Query/View permission deny - strange behaviour

[Birgit Laggner]

Hi Pietro,

you try to grant permissions to user catasto_admins - maybe the problem 
lies here? You don't write about permission settings for catasto_admins, 
so I can't be sure.


Regards,

Birgit

Am 11.05.2016 um 12:39 schrieb Pietro Rossin:

Hello
I have problems with postgres/postgis query/view permissions.

Postgresql 9.4 Postgis 2.1

I created a user pippo inherit.

I created a group catasto_editors inherit and a group basi_viewers.

I granted catasto_editors to pippo and basi_viewers to catasto_editors.

Then two schemas basi and catasto

GRANT USAGE ON SCHEMA catasto TO catasto_editors;
GRANT USAGE ON SCHEMA basi TO basi_viewers;

All tables/views in catasto are selectable/editable/updateable ecc ecc by
catasto_editors.
All tables/views in basi are selectable by basi_viewers.


Then I made a view in catasto

***
CREATE OR REPLACE VIEW catasto.points_in_polygons AS
  SELECT a.id,
 b.name,
 b.code,
 a.id_impianto,
 a...
 FROM catasto.points a,
 basi.polygons b
   WHERE st_contains(b.geom, a.geom);

ALTER TABLE catasto.points_in_polygons   OWNER TO catasto_admins;
GRANT ALL ON TABLE catasto.points_in_polygons TO catasto_admins;
GRANT SELECT, UPDATE, INSERT, DELETE, TRIGGER ON TABLE
catasto.points_in_polygons TO catasto_editors;


This doesn't work, an error occurs saying "ERROR: permission denied for
relation polygons SQL state: 42501"

Strange thing, if I execute:

**
SELECT a.id,
 b.name,
 b.code,
 a.id_impianto,
 a...
 FROM catasto.points a,
 basi.polygons b
   WHERE st_contains(b.geom, a.geom);
***

the query itself works...

Also if I execute
*
SELECT count(id)
   FROM basi.polygons;
*

I get the item count..

Where is my mistake?

Thanks
Pietro









--
View this message in context: 
http://postgis.17.x6.nabble.com/Query-View-permission-deny-strange-behaviour-tp5010068.html
Sent from the PostGIS - User mailing list archive at Nabble.com.
___
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] Average Algorithm for Raster Resampling

[Birgit Laggner]

Dear list,

is there anybody who can tell me if the algorithm Average is available 
for Raster Resampling? It is not listed as an algorithm in the raster 
reference of the PostGIS manual but somehow I am hoping it might still 
be implemented :-) We are using PostGIS version 2.2.1 
(postgis_full_version(): "POSTGIS="2.2.1 r14555" GEOS="3.5.0-CAPI-1.9.0 
r4084" PROJ="Rel. 4.9.1, 04 March 2015" GDAL="GDAL 1.11.2, released 
2015/02/10" LIBXML="2.9.2" LIBJSON="0.11.99" TOPOLOGY RASTER").


Thanks and regards,

Birgit


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

Re: [postgis-users] A function for “Esri union” union on big tables on github.

[Birgit Laggner]

Hi Lars,

awesome, thank you. I will surely test your code at the next occasion.

Regards,

Birgit

Am 10.02.2016 um 09:11 schrieb Lars Aksel Opsahl:

Hi


There has been different mails about this topic lately. We have now have added 
the code we use to Github and hopefully somebody can pick up some ideas or just 
use this function as it is.


The basic idea is that you call this function and with 2 tables as input. The 
following happens in the function

   *   Builds up a content based grid

   *   Computes the result

   *   Removes the grid lines from the result

   *   Returns a table name with the union of this two tables. For areas that 
intersect you get attributes from both tables and for areas that only exits in 
one of the tables you only get attributes from one table.


The code is found at https://github.com/larsop/esri_union


About performance. The code added now runs in a single thread, but we have a 
slightly modified code that runs in parallel using “Gnu parallel” and then we 
can increase the performance many times depending on how many CPU you have on 
your server. Here is an example running with 20 threads.


num points num polygons table size

Table 1 40435700 1088614 637 MB

Table 1 933145431 7924019 10127 MB

Result table 2042294001 43668256 30 GB


The time used to do the intersection was 152 minutes. I will add the parallel 
code later when I have time to make the code ready.


Lars

___
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] Handling LINESTRINGM

[Birgit Laggner]

Hi Paolo,

perhaps converting into ewkt or ewkb format with st_asewkt() or 
st_asewkb() could be a workaround? I haven't worked with geometries 
including m-coordinates yet, so this is only an idea...


Good luck and regards,

Birgit

Am 22.12.2015 um 13:27 schrieb Paolo Cavallini:

Hi all,
I discovered that LINESTRINGM are lost during a dump/restore (the table
is created, but completely empty). Is this a known bug?
So I tried converting them to LINESTRING, but could not find a
reasonable way: anyone has a suggestion?
Thanks.


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

Re: [postgis-users] Errors restoring a db

[Birgit Laggner]

You don't have to use postgres user per se - but the user using the 
'copy' command has to be a superuser. Perhaps it might be easier to 
change your own user role accordingly?


Am 22.12.2015 um 14:25 schrieb Paolo Cavallini:

Il 22/12/2015 14:21, Humberto Cereser Ibanez ha scritto:

Hi Paolo,

If your sql is trying to use 'copy' command, I experencied that the
superuser postgres can run it and stop the errors that you showed.

I see, thanks. Not very nice to be forced using postgres though.
All the best.


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

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 

*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 >
*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 

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] Problem updating geometry in a record

[Birgit Laggner]

Hi Bob,

I had the problem sometimes that st_union() result was NULL. Perhaps 
that's the case with your query, too? Is it necessary to use st_union()? 
Otherwise you could replace it with st_collect(), if you just want to 
collect all geometries into a multilinestring. That's more robust 
because no geoprocessing is needed...


Regards,

Birgit


Am 28.09.2015 um 22:46 schrieb Bistrais, Bob:


I have been able to derive a geometry from a query.  I want to insert 
that geometry into a Multilinestring field.  I have some SQL code to 
do it.  When I run the code, the output say that the query was 
successful. But I don’t see the geometry inserted into the table.  
Here is the code:


update highschools_aroo set routegeom =

(SELECT ST_Union(geom) as mygeom from (SELECT pt.geom

FROM pgr_dijkstra(

'SELECT gid AS id,

sourceroute::integer as source,

targetroute::integer as target,

(feet / 5280)::double precision AS cost,

feet::double precision as reverse_cost

FROM ngrdsaroo2',2465,1713,false,true)

as di

JOIN ngrdsaroo2 pt on di.id2 = pt.gid) as foo)

WHERE source_node = 2465 and target_node = 1713;

-As mentioned, this runs with no errors, but the geometry field is not 
updated.  Any suggestions?



Thanks,

Bob



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


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

Re: [postgis-users] Trouble with returning MultiPolygon on psql query

[Birgit Laggner]

Hi Joseph,

I tried your intersects query after creating the multipolygon using the 
st_geomfromtext function with the coordinates from your geojson objects. 
There were no problems with the intersects query. The query returned a 
result even with the point within the lower polygon.


So my question would be: Are you sure, you did import all multipolygon 
parts in postgres? If you can't visualize the PostGIS geometries, 
perhaps you could dump the multipolygon in question to single polygons 
and compare the resulting number with the original as a first approach?


Regards,

Birgit


Am 12.09.2015 um 19:29 schrieb Joseph Spenner:
I believe I may have overcomplicated my question.  My apologies.  I've 
come up with a much simpler explanation of the issue:


I have a MultiPolygon with 2 relatively simple polygons in it:

http://microflush.org/json/MultiPolygon.json 




I've pulled out the 2 polygons from the above MultiPolygon below:

http://microflush.org/json/upper.json

http://microflush.org/json/lower.json

Any/all of the 3 above can be copied/pasted into a GeoJSON tester to 
view them:

http://geojsonlint.com/


I've stored the original MultiPolygon in PostgreSQL as described in my 
original post:


$ ogr2ogr -f "PostgreSQL" PG:"dbname=weatherzones user=postgres" 
"MultiPolygon.json" -nln polys


Here's how it looks in Postgres after the import:

polytest=# \d polys;
   Table "public.polys"
Column|  Type | Modifiers
--+-+-
 ogc_fid  | integer | not null default 
nextval('polys_ogc_fid_seq'::regclass)

 wkb_geometry | geometry(Geometry,4326) |
 warnings | character varying   |
Indexes:
"polys_pk" PRIMARY KEY, btree (ogc_fid)
"polys_geom_idx" gist (wkb_geometry)

polytest=#


When I try to query postgres using single points which lie in the 
lower polygon, I do not get rows returned.
However, when I query using points which lie in the upper polygon, I 
get the row returned.


ie:

This point lies in the lower polygon, and this query returns no rows:

polytest=# select warnings from polys where 
ST_Intersects(ST_PointFromText('POINT( -116.024551 38.485773 )', 
4326), wkb_geometry);



This next point lies in the upper polygon, and this query returns a 
row, which is the MultiPolygon:


polytest=# select warnings from polys where 
ST_Intersects(ST_PointFromText('POINT( -114.879913 39.249129 )', 
4326), wkb_geometry);


Is there something wrong with my query, or perhaps the original 
MultiPolygon which can explain why a query with a point from the lower 
polygon won't return the MultiPolygon row?  I've done this with other 
MultiPolygons and have not had an issue like this.


Any help would be great.

Thanks!

Joseph Spenner






*From:* Joseph Spenner 
*To:* PostGIS Users Discussion 
*Sent:* Wednesday, September 9, 2015 12:10 PM
*Subject:* Trouble with returning MultiPolygon on psql query

Hello, I have a file (source.json) containing 25 GeoJSON polygons with 
weather alerts, which I imported into postgres using the following 
command:


$ ogr2ogr -f "PostgreSQL" PG:"dbname=weatherzones user=postgres" 
"source.json" -nln polys


Here's how it looks in Postgres:

polytest=# \d polys;
Table "public.polys"
Column|  Type | Modifiers
--+-+-
 ogc_fid  | integer | not null default 
nextval('polys_ogc_fid_seq'::regclass)

 wkb_geometry | geometry(Geometry,4326) |
 warnings | character varying |
Indexes:
"polys_pk" PRIMARY KEY, btree (ogc_fid)
"polys_geom_idx" gist (wkb_geometry)

polytest=#

I can see it stored all of them:

polytest=# select count(*) from polys;
 count
---
25
(1 row)

polytest=#

However, when I try to use a ST_Intersects(ST_PointFromText('POINT(), 
I can't seem to get lines with MultiPolygons to return rows.  The 
regular Polygons seem to work fine, though.


For example, I have the following MultiPolygon in my table from Nevada:
http://microflush.org/json/24.json

Using a validator site, such as http://geojsonlint.com/ , I can 
copy/paste that JSON into the interface and it looks like 2 adjacent 
shapes making up the MultiPolygon.


If I grab points from the upper section of the polygon, using 
something like http://itouchmap.com/latlong.html, I get rows returned:
polytest=# select warnings from polys where 
ST_Intersects(ST_PointFromText('POINT( -114.543457 39.336281 )', 
4326), wkb_geometry) limit 2;
polytest=# select warnings from polys where 
ST_Intersects(ST_PointFromText('POINT( -115.658569 39.500297 )', 
4326), wkb_geometry) limit 2;


But if I use a point from the lower Polygon section of that 
MultiPolygon, I get zero:

Re: [postgis-users] Question on optimizing slow geospatial query

[Birgit Laggner]

Hi Trang,

I think, it could help to create btree indices on startts and uuid, 
too, since you are using them in your where clause as a filter (a 
probably unnecessary question regarding your date filter: I would expect 
the result of t.startts  '2015-01-16' and t.startts  '2015-01-17' to 
be null because there are no days between the 16th and the 17th of 
january - but perhaps it was only an example...). And in general, my 
suggestion would be to reduce the use of st_intersects to the necessary 
mimimum.
You could use the with clause for the filtration of your input data and 
afterwards double join the two tables first on the startloc and second 
on the endloc for the assignment of the origin and the destination zone. 
Then group by origin and destination zones while counting your trips and 
you should have your end result.


Here is how I would imagine the query:

with
t as (select trip_id, startloc, endloc from od1.trip_v1 where startts 
between 'minimum start date' and 'maximum start date'),
z as (select zone from od1.taz where uuid in ('kansas_303', 
'kansas_601', 'kansas_603', etc))


select  count(t.trip_id) as number_of_trips, orig.zone as orig_zone, 
dest.zone as dest_zone from t left join z as orig on 
st_intersects(t.startloc, z.geom) left join z as dest on 
st_intersects(t.endloc, z.geom) group by orig.zone, dest.zone;



However, I am not sure about how the gist indices work together with the 
subselects of the with clause...


Hope this helps,

Birgit.



Am 11.02.2015 um 08:19 schrieb Trang Nguyen:


Hi,

I am a newbie to Postgres/PostGIS and have a long running query that I 
would like to optimize.


There are two tables (trip and zone) that I am joining in the query, 
one which has “startloc” and “endloc” columns with type 
Geometry(Point) and other which contains a Geometry(MultiPolygon). 
There are GIST indexes on all above columns:


CREATE TABLE od1.trip_v1

(

  pkey bigint NOT NULL,

  trip_id character varying,

  startts timestamp without time zone,

  endts timestamp without time zone,

  startloc geometry(Point),

  endloc geometry(Point),

  …

)

CREATE TABLE od1.taz

(

  uuid character varying NOT NULL,

  zone character varying NOT NULL,

  createdts timestamp without time zone NOT NULL,

  updatedts timestamp without time zone NOT NULL,

  geom geometry(MultiPolygon) NOT NULL,

  CONSTRAINT taz_pkey PRIMARY KEY (uuid)

)


I’m interested in building a matrix that, for a given set of input 
zones, returns trips along with their start and end zones. Output 
looks like:


10 trips that start at Zone A, ends at Zone B

2 trips that start at Zone A, ends at Zone C

9 trips that start at Zone A, ends at other

13 trips that start at Zone C, ends at Zone D

Since I am dealing with a large dataset ( 24 million records and 
growing), I was planning on writing a query that returns the trips 
grouped by each zone along with match condition (start, end or both) 
and doing aggregation on the client layer. I’m not sure whether this 
is the best approach but I expect that otherwise, I would end up 
having to write a very complex query to handle that type of aggregation.


Even so, the current query is very slow with very high cost:

SELECT t.trip_id,

case

when ST_intersects(t.startloc, z.geom) and 
ST_intersects(t.endloc, z.geom) then 'orig-dest'


when ST_intersects(t.startloc, z.geom) then 'orig'

when ST_intersects(t.endloc, z.geom) then 'dest'

else 'none'

end  as match_cond,

z.zone from od1.trip_v1   t, od1.taz z

where t.startts  '2015-01-16' and t.startts  '2015-01-17'

and z.uuid in ('kansas_303', 'kansas_601', 'kansas_603', etc)

and ST_intersects(t.startloc, z.geom)

or ST_intersects(t.endloc, z.geom)

group by z.zone, t.trip_id, match_cond;

Explain plan:

Group (cost=231446695055.73..245971533247.59 rows=14240037443 
width=3498)


  -  Sort (cost=231446695055.73..231482295149.34 rows=14240037443 
width=3498)


Sort Key: z.zone, t.trip_id, (CASE WHEN ((t.startloc  
z.geom) AND _st_intersects(t.startloc, z.geom) AND (t.endloc  
z.geom) AND _st_intersects(t.endloc, z.geom)) THEN 'orig-dest'::text 
WHEN ((t.startloc  z.geom) AND _st_intersects(t.startloc, (...)


-  Nested Loop (cost=91.70..14401634128.24 rows=14240037443 
width=3498)


  -  Seq Scan on taz z (cost=0.00..739.19 rows=4619 
width=3406)


  -  Bitmap Heap Scan on trip_v1 t  (cost=91.70..4151.26 
rows=453 width=107)


Recheck Cond: ((startloc  z.geom) OR (endloc  
z.geom))


Filter: (((startts  '2015-01-16 
00:00:00'::timestamp without time zone) AND (startts  '2015-01-17 
00:00:00'::timestamp without time zone) AND ((z.uuid)::text = ANY 
('{kansas_303,kansas_601,kansas_603,kansas_604,kansas_10,kansas_11,kan 
(...)


-  BitmapOr (cost=91.70..91.70 rows=2706 width=0)

  -  Bitmap Index Scan on 
idx_trip_v1_startloc  

Re: [postgis-users] PostGIS supported by which ArcGIS

[Birgit Laggner]

Hi David,

I have no experiences with ArcSDE, but I can say, that no dba had to 
esrify our database in order to enable the interoperability between 
ArcGIS and PostgreSQL/PostGIS. We have a simple PostgreSQL-DB (version 
9.1) with an PostGIS extension (version 2.1). All we had to do was 
download a bunch of files from the ArcGIS Helpdesk, copy them into the 
appropriate ArcGIS program folder of our locally installed desktop 
ArcGIS, and then everything worked. We had to define the connection to 
the PostgreSQL database in ArcGIS and then we were able to access the 
tables of our PostgreSQL database including the PostGIS geometries (not 
the raster data!) via read/write access.


If ArcGIS writes geometries into the PostgreSQL database it names the 
geometry column as shape with the datatype geometry - no pg_geometry 
or st_geometry.


Regards,

Birgit.


Am 08.10.2014 17:28, schrieb David Fawcett:
I assume that on the server side, your dba has 'Esrified' the database 
by using a one of the ArcGIS desktop tools to convert your existing 
Postgres database into an 'ESRI Enterprise Geodatabase, which creates 
the sde schema and all of the sde functions, etc.  This costs an 
ArcGIS Server license.


We are doing this for our primary production spatial database.  One 
thing that we are doing is using the PG_GEOMETRY option to specify the 
PostGIS spatial datatype instead of ESRI's ST_GEOMETRY datatype.  This 
allows us to use all of the PostGIS functionality on features created 
with ESRI tools.  All of my backend processes are written in 
Python/SQL/PgPSQL without any ESRI tools.


I would be very surprised if you can edit geometries in PostGIS using 
the ArcGIS tools without using the SDE functionality and licensing.  
(I would love it if you can prove me wrong.)  Note that the SDE 
functionality is now built into all of ESRI's products, and the SDE 
server software will be depricated after 10.2.




On Thu, Oct 2, 2014 at 2:30 AM, Birgit Laggner 
birgit.lagg...@ti.bund.de mailto:birgit.lagg...@ti.bund.de wrote:


We use ArcGIS 10.2 for read/write interaction with our PostgreSQL
database. We only have ArcGIS Desktop without SDE. There is a
driver (database client file) from ArcGIS for PostgreSQL, you have
to install. We are more or less content with this concept. Surely,
the performance is not as if you would be working directly on the
PostgreSQL database, but still, it isn't necessary to always
export and import into and from shapefile. That's a huge plus in
my opinion.

Here is a link to a documentation site of ArcGIS regarding the
requirements for the interoperability with PostgreSQL:


http://resources.arcgis.com/en/help/system-requirements/10.2/index.html#//0151007500

Regards,

Birgit.


Am 27.09.2014 20 tel:27.09.2014%2020:28, schrieb Stefan Keller:

Actually yes, I meant the PostGIS option of storing data in
PostgreSQL/PostGIS from ArcGIS.

--S.

2014-09-27 20:20 GMT+02:00 Randal Hale
rjh...@northrivergeographic.com
mailto:rjh...@northrivergeographic.com:

as far as I understand it (and I could be wrong on some of
this)

Desktop will support reading if you load the postgresql
library files (that
can be downloaded form your ESRI customer care portal).
Once those are loaded you can read through a definition
query data from
postgis into arcgis. I just made a map using arcgis and
reading data from
postgis.
It doesn't support write unless you have SDE for
Postgresql - which isn't
postgis. That dives into the world of ArcGIS Server


ArcGIS supposedly is moving to being able to read/write
the geopackage
format which to me opens read/write to PostGIS. BUT if
they open read write
for postgis that starts a slow death of SDEso I don't
believe it ever
will (for the foreseeable future).

Randy




On 09/27/2014 02:10 PM, Stefan Keller wrote:

I'm asked from time to time for an advice which ArcGIS product
supports PostGIS (read/write)?
Even after reading the Esri pages I'm not sure if it's
only ArcGIS
for Server Workgroup (formerly ArcSDE).
Does it work also with ArcGIS for Desktop Basic (ArcView)?

Yours, S.
___
postgis-users mailing list
postgis-users@lists.osgeo.org
mailto:postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users


--
-
Randal Hale
North River Geographic Systems, Inc
http://www.northrivergeographic.com

Re: [postgis-users] PostGIS supported by which ArcGIS

[Birgit Laggner]

We use ArcGIS 10.2 for read/write interaction with our PostgreSQL 
database. We only have ArcGIS Desktop without SDE. There is a driver 
(database client file) from ArcGIS for PostgreSQL, you have to install. 
We are more or less content with this concept. Surely, the performance 
is not as if you would be working directly on the PostgreSQL database, 
but still, it isn't necessary to always export and import into and from 
shapefile. That's a huge plus in my opinion.


Here is a link to a documentation site of ArcGIS regarding the 
requirements for the interoperability with PostgreSQL:


http://resources.arcgis.com/en/help/system-requirements/10.2/index.html#//0151007500

Regards,

Birgit.


Am 27.09.2014 20:28, schrieb Stefan Keller:

Actually yes, I meant the PostGIS option of storing data in
PostgreSQL/PostGIS from ArcGIS.

--S.

2014-09-27 20:20 GMT+02:00 Randal Hale rjh...@northrivergeographic.com:

as far as I understand it (and I could be wrong on some of this)

Desktop will support reading if you load the postgresql library files (that
can be downloaded form your ESRI customer care portal).
Once those are loaded you can read through a definition query data from
postgis into arcgis. I just made a map using arcgis and reading data from
postgis.
It doesn't support write unless you have SDE for Postgresql - which isn't
postgis. That dives into the world of ArcGIS Server


ArcGIS supposedly is moving to being able to read/write the geopackage
format which to me opens read/write to PostGIS. BUT if they open read write
for postgis that starts a slow death of SDEso I don't believe it ever
will (for the foreseeable future).

Randy




On 09/27/2014 02:10 PM, Stefan Keller wrote:

I'm asked from time to time for an advice which ArcGIS product
supports PostGIS (read/write)?
Even after reading the Esri pages I'm not sure if it's only ArcGIS
for Server Workgroup (formerly ArcSDE).
Does it work also with ArcGIS for Desktop Basic (ArcView)?

Yours, S.
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users


--
-
Randal Hale
North River Geographic Systems, Inc
http://www.northrivergeographic.com
423.653.3611 rjh...@northrivergeographic.com
twitter:rjhale http://about.me/rjhale
http://www.northrivergeographic.com/spatial-connect

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



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

Re: [postgis-users] Spike finder

[Birgit Laggner]

Hi L.

I would suspect the first assignment of the result_pnt (result_pnt := 
ST_Collect(result_pnt, st_pointn(lineusp, point_id));) as the cause of 
the problem: Within your assignment, you try to use the geometry of 
result_pnt which does not yet exist at this moment. Maybe, you could 
work with an IF clause to assign some sort of start value to your 
result_pnt variable.


Regards,

Birgit.


Am 11.06.2014 10:56, schrieb Pier Lorenzo Marasco:
I tried to modify spikeRemover from Schmidt  Krüger to obtain just 
the spike position. For some reason, that I can't figure out, my 
script fails and I don't obtain a multipoint position. Probably there 
is a problem in the aggregation of different points but I don't know 
how to manage it in another way. I tried with ST_collect and ST_union 
but with both I don't obtain a right geometry.


CREATE OR REPLACE FUNCTION ST_SpikeFinder (geometry, angle double 
precision)

 returns geometry as
 $body$
 DECLARE
 ingeom alias for $1;
 angle  alias for $2;
 lineusp geometry;
 newgeom geometry;
 numpoints integer;
 point_id integer;
 result_pnt geometry;

begin
-- input geometry or rather set as default for the output
newgeom := ingeom;
   -- check polygon
if (select st_geometrytype(ingeom)) = 'ST_Polygon' then
if (select st_numinteriorrings(ingeom)) = 0 then
lineusp := st_boundary(ingeom) as line;
numpoints := st_numpoints(lineusp);
point_id := 0;
-- the geometry passes pointwisely
while (point_id = numpoints) loop
result_pnt := ST_Collect(result_pnt, 
st_pointn(lineusp, point_id));
-- the check of the angle at the current point of 
a spike including the special case, that it is the first point.
if (select abs(pi() - 
abs(st_azimuth(st_pointn(lineusp, case when point_id= 1 then 
st_numpoints(lineusp) - 1 else point_id - 1 end),
st_pointn(lineusp, point_id)) - 
st_azimuth(st_pointn(lineusp, point_id), st_pointn(lineusp, point_id + 
1) = angle then
(probably the problem is here)  ---  result_pnt := 
ST_Union(result_pnt, st_pointn(lineusp, point_id));

end if;
point_id = point_id + 1;
end loop;
end if;
end if;
return result_pnt;
 end;
 $body$
   language 'plpgsql' volatile;

Any idea ?
Thank you in advance.


L.


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


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

Re: [postgis-users] ST_DWithin advice

[Birgit Laggner]

Hi Adam,

perhaps you could solve the problem with this workaround:

1. Transform pt1 into geometry data type with srid 4326 (coordinates are 
in degrees, now)

2. Create a new point by shifting pt1 for # degrees in x-direction
3. Create a linestring geometry by connecting pt1 and the new point
4. Transform the linestring back into geography data type
5. Calculate length of linestring in meters
6. Use length in meters for definition of the ST_DWithin radius

It's a little bit complicated and maybe you came up with a better idea 
in the meantime...


Regards,

Birgit.



Am 02.06.2014 20:58, schrieb Adam Wright:
We have several tables with latitude and longitude columns (data type: 
numeric) and I need to calculate whether a given lat/lon is within # 
degrees of a point-radius ring. I came up with the function below. The 
user supplied radius has to be in degrees (e.g. give me all records 
within 5 to 25 degrees of this lat/lon pair). Sample query:  select * 
from mytable where 
dist_check(mytable.latitude.mytable.longitude,35.6895,139.6917,5,25)=1.


Any advice on solving the same problem using the geography data type 
when the input radius is supplied in degrees?


CREATE OR REPLACE FUNCTION dist_check(lat1 numeric, lon1 numeric, lat2 
numeric, lon2 numeric, innerradius numeric, outerradius numeric)


pt1 geometry;
pt2 geometry;
BEGIN
pt1 := ST_MakePoint(lon1,lat1);
pt2 := ST_MakePoint (lon2,lat2);
IF ST_DWithin(pt1,pt2,outerRadius) AND NOT ST_DWithin(pt1,pt2,innerRadius)
THEN return 1;
ELSE
return 0;
END IF;

Thanks!


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


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

Re: [postgis-users] FYI: Solution to Using PostGIS in ArcGIS.

[Birgit Laggner]

Hi,

we work with ArcGIS 10.1, now, and since, we don't need ArcSDE for 
communication with PostgreSQL db and PostGIS geometries. If you install 
special db client libraries for PostgreSQL, you are able to create a 
direct database connection to the PostgreSQL database and to directly 
im- and export data between ArcGIS and PostgreSQL. We have not much 
improved our testing because we are waiting for SP1 for ArcGIS 10.1 
(which should support PostgreSQL 9.1.3 and PostGIS 2.0) and the 
installation of our new PostgreSQL server, so e.g. I can't say which 
ArcGIS licence level you would need...


For further information see:

PostgreSQL database requirements for ArcGIS 10.1:
http://resources.arcgis.com/en/help/system-requirements/10.1/index.html#//0151007500

Setting up a connection to PostgreSQL:
http://resources.arcgis.com/en/help/main/10.1/index.html#//002p003q00

Regards,

Birgit.


Am 07.11.2012 09:25, schrieb Paolo Corti:

On Wed, Nov 7, 2012 at 7:39 AM, Robert Buckleyrobertdbuck...@yahoo.com  wrote:

Hi,

just my 2Cents..ArcGIS is now realising that there is a great need for
FOSSGIS interfaces within the ESRI Architecture. If you have ArcSDE and
ArcEditor (only a few 10s of thousands) you can work with PostgreSQL and
PostGIS dbs.
Other than that, there may be some opensource plugins or extentions that you
may be able to test.

http://gis.stackexchange.com/questions/147/how-can-i-connect-to-a-postgis-database-from-arcmap-9-3-and-10-0
http://www.paolocorti.net/2008/06/06/spatial-database-for-postgres-and-arcgis-users-how-to-choose/


Hello Robert

please be aware that zigGIS was discontinued (the posted links are
very old) [1].
In any case, the source code has always been around [2] and I think it
should be quiete easy to adapt it to work to latest ArcGIS versions.

regards
p

[1] https://groups.google.com/forum/?fromgroups=#!topic/ziggis/FK6KNRVk3p4
[2] http://code.google.com/p/ziggis/

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users


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