[postgis-users] Finding closet intersect point along direction of line

[adrian.ki...@dse.vic.gov.au]

Hi All

I'm looking to see if anyone has a PostGIS solution to finding the closest point on another polyline to a line end point along the line's azimuth. The picture attached shows what I'm looking for. For end point A the usual closest point is B but I want C. The picture shows a polygon but I expect a conversion to polyline to get points on the edge.

Cheers

Adrian

Using: PostGIS 2.0 on PostgreSQL 9.2, Windows 7

Notice:
This email and any attachments may contain information that is personal, confidential,
legally privileged and/or copyright. No part of it should be reproduced, adapted or communicated without the prior written consent of the copyright owner.

It is the responsibility of the recipient to check for and remove viruses.

If you have received this email in error, please notify the sender by return email, delete it from your system and destroy any copies. You are not authorised to use, communicate or rely on the information contained in this email.

Please consider the environment before printing this email.

[Simon Greener]

Adrian,

Below is a function that I found on the web that should do what you want.

CREATE OR REPLACE FUNCTION ST_Extend(eje_ geometry, bound_ geometry)
RETURNS geometry
AS $$
-- Use Case: I need to "extend" a linestring in another one. The new linestring must be
-- obtained such as their new extremes are the intersection of it with a
-- polygon that contains it. The original linestring is composed of only a
-- linear segment.
-- version: alfa , by Julio A. Galindo, April 17, 2007: juliogalindoq at gmail.com
DECLARE
b_ geometry = boundary(bound_);
dist float;
max_dist float = 0;
n_points int;
pto_1 geometry;
pto_2 geometry;
first_pto geometry;
last_pto geometry;
u_1 float;
u_2 float;
norm float;
result text = 'LINESTRING(';
BEGIN
IF ( GeometryType(eje_) NOT LIKE 'LINESTRING' OR
GeometryType(bound_) NOT LIKE 'POLYGON' THEN
RETURN NULL;
END IF;
-- First Search how far is the boundary: (worst case)
pto_1 := StartPoint(eje_);
pto_2 := EndPoint(eje_);
FOR i IN 1..NumPoints(b_)-1 LOOP
dist := distance(PointN(b_,i),pto_1);
IF dist > max_dist THEN max_dist := dist; END IF;
max_dist := dist;
END IF;
dist := distance(PointN(b_,i),pto_2);
IF dist > max_dist THEN max_dist := dist; END IF;
max_dist := dist;
END IF;
END LOOP;
-- Now extent the linestring:
pto_2 := PointN(eje_,2);
u_1 := X(pto_2)-X(pto_1);
u_2 := Y(pto_2)-Y(pto_1);
norm := sqrt(u_1^2 + u_2^2);
first_pto := MakePoint(X(pto_1)-u_1/norm*dist,Y(pto_1)-u_2/norm*dist);
n_points := nPoints(eje_);
IF n_points > 2 THEN
pto_1 := PointN(eje_,n_points-1);
pto_2 := PointN(eje_,n_points);
u_1 := X(pto_2)-X(pto_1);
u_2 := Y(pto_2)-Y(pto_1);
norm := sqrt(u_1^2 + u_2^2);
END IF;
last_pto := MakePoint(X(pto_2)+u_1/norm*dist,Y(pto_2)+u_2/norm*dist);
result := result || X(first_pto) || ' ' || Y(first_pto) || ',';
FOR i IN 1..NumPoints(eje_) LOOP
result := result || X(PointN(eje_,i)) || ' ' || Y(PointN(eje_,i)) || ',';
END LOOP;
result := result || X(last_pto) || ' ' || Y(last_pto) || ')';
-- Find the final Linestring:
b_ := intersection(GeomFromText(result,SRID(eje_)),bound_);
RETURN b_;
END $$
LANGUAGE plpgsql
STABLE
RETURNS NULL ON NULL INPUT;

-- Test 1:
select asText(extendLine(geomFromText('LINESTRING(5 5, 12 12)',-1),GeomFromText('POLYGON((0 0,15 0,15 20,0 20,0 0))',-1)));

select asText(extendLine(geomFromText('LINESTRING(5 5, 12 12, 13 12)',-1),GeomFromText('POLYGON((0 0,15 0,15 20,0 20,0 0))',-1)));

regards
Simon

On Wed, 10 Apr 2013 15:34:04 +1000, <adrian.ki...@dse.vic.gov.au> wrote:

> Hi All
>
> I'm looking to see if anyone has a PostGIS
> solution to finding the closest point on another polyline to a line end
> point along the line's azimuth. The picture attached shows what I'm looking
> for. For end point A the usual closest point is B but I want C. The picture
> shows a polygon but I expect a conversion to polyline to get points on
> the edge.
>
> Cheers
>
> Adrian
>
> Using: PostGIS 2.0 on PostgreSQL 9.2,
> Windows 7
>
>
>
>
>
> Notice:

> This email and any attachments may contain information that ispersonal, confidential,
> legally privileged and/or copyright. No part of itshould be reproduced, adapted or communicated without the prior written consentof the copyright owner.

>
>
> It is the responsibility of the recipient to check for and removeviruses.
>
>
>

> If you have received this email in error, please notify the sender by returnemail, delete it from your system and destroy any copies. You are not authorisedto use, communicate or rely on the information contained in this email.

>
>
>
> Please consider the environment before printing this email.
>
>

--
Holder of "2011 Oracle Spatial Excellence Award for Education and Research."
SpatialDB Advice and Design, Solutions Architecture and Programming,
Oracle Database 10g Administrator Certified Associate; Oracle Database 10g SQL Certified Professional
Oracle Spatial, SQL Server, PostGIS, MySQL, ArcSDE, Manifold GIS, FME, Radius Topology and Studio Specialist.
39 Cliff View Drive, Allens Rivulet, 7150, Tasmania, Australia.
Website: www.spatialdbadvisor.com
Email: si...@spatialdbadvisor.com
Voice: +61 362 396397
Mobile: +61 418 396391
Skype: sggreener
Longitude: 147.20515 (147° 12' 18" E)
Latitude: -43.01530 (43° 00' 55" S)
GeoHash: r22em9r98wg
NAC:W80CK 7SWP3
_______________________________________________
postgis-users mailing list
postgi...@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

[adrian.ki...@dse.vic.gov.au]

Thanks for the code Simon. After a few tweaks I got it working perfectly in PostGIS2.

Here's the updated code.

Adrian

CREATE OR REPLACE FUNCTION ST_Extend(eje_ geometry, bound_ geometry)
RETURNS geometry
AS $$
-- Use Case: I need to "extend" a linestring in another one. The new linestring must be
-- obtained such as their new extremes are the intersection of it with a
-- polygon that contains it. The original linestring is composed of only a
-- linear segment.
-- version: alfa , by Julio A. Galindo, April 17, 2007: juliogalindoq at gmail.com
DECLARE

      b_ geometry = st_boundary(bound_);

      dist float;
      max_dist float = 0;
      n_points int;
      pto_1 geometry;
      pto_2 geometry;
      first_pto geometry;
      last_pto geometry;
      u_1 float;
      u_2 float;
      norm float;
      result text = 'LINESTRING(';
BEGIN

      IF  st_GeometryType(eje_)   NOT LIKE 'ST_LineString' OR
           st_GeometryType(bound_) NOT LIKE 'ST_Polygon' THEN

         RETURN NULL;
      END IF;
      -- First Search how far is the boundary: (worst case)

      pto_1 := st_StartPoint(eje_);
      pto_2 := st_EndPoint(eje_);
      FOR i IN 1..st_NumPoints(b_)-1 LOOP
         dist := st_distance(st_PointN(b_,i),pto_1);
         IF dist > max_dist THEN max_dist := dist; --END IF;

            max_dist := dist;
         END IF;

         dist := st_distance(st_PointN(b_,i),pto_2);
         IF dist > max_dist THEN max_dist := dist; --END IF;

            max_dist := dist;
         END IF;
      END LOOP;
      -- Now extent the linestring:

      pto_2 := st_PointN(eje_,2);
      u_1 := st_X(pto_2)-st_X(pto_1);
      u_2 := st_Y(pto_2)-st_Y(pto_1);

      norm := sqrt(u_1^2 + u_2^2);

      first_pto := st_MakePoint(st_X(pto_1)-u_1/norm*dist,st_Y(pto_1)-u_2/norm*dist);
      n_points := st_nPoints(eje_);
      IF n_points > 2 THEN
         pto_1 := st_PointN(eje_,n_points-1);
         pto_2 := st_PointN(eje_,n_points);
         u_1   := st_X(pto_2)-st_X(pto_1);
         u_2   := st_Y(pto_2)-st_Y(pto_1);

         norm  := sqrt(u_1^2 + u_2^2);
      END IF;

      last_pto := st_MakePoint(st_X(pto_2)+u_1/norm*dist,st_Y(pto_2)+u_2/norm*dist);
      result   := result || st_X(first_pto) || ' ' || st_Y(first_pto) || ',';
      FOR i IN 1..st_NumPoints(eje_) LOOP
          result := result || st_X(st_PointN(eje_,i)) || ' ' || st_Y(st_PointN(eje_,i)) || ',';
      END LOOP;
      result := result || st_X(last_pto) || ' ' || st_Y(last_pto) || ')';

      -- Find the final Linestring:

      b_ := st_intersection(st_GeomFromText(result,st_SRID(eje_)),bound_);

      RETURN b_;
  END $$
  LANGUAGE plpgsql
  STABLE
  RETURNS NULL ON NULL INPUT;

Notice:
This email and any attachments may contain information that is personal, confidential,
legally privileged and/or copyright. No part of it should be reproduced, adapted or communicated without the prior written consent of the copyright owner.

It is the responsibility of the recipient to check for and remove viruses.

If you have received this email in error, please notify the sender by return email, delete it from your system and destroy any copies. You are not authorised to use, communicate or rely on the information contained in this email.

[Simon Greener]

Adrian,

That's excellent.

regards

Simon

[pcr...@pcreso.com]

Hi, I'm trying to determine all the possible intersections in sets of polygons. The approach I'm using is to load each set of polygon perimeters into a Postgis topology & extract all the faces as geometries. This results in some new polygons lying between the input ones. Which is the nature of this technique- that is not the problem. I then determine the number of original polygons that the new ones overlap, which should go from 0 up. My problem arises in that I'm getting new polygons, derived from the boundaries of the old ones, which are actually between old polygons, but have a non-zero (but very small - like 1e-15) area of overlap - but all they have in common is the boundary line which should have a zero area. So some of my new polygons supposedly overlap an old polygon, but in fact they do (or should) not. My solution at the moment is to reject (treat as a zero area) any overlaps < 1e-15 in area, so the number of overlaps represents the "real" data. It works for my purposes, but is not ideal. Is this a bug? An inevitable result with near zero values & finite precision? Is there a better way of dealing with this issue?  Thanks,     Brent Wood

[Sandro Santilli]

On Tue, Apr 16, 2013 at 05:05:09PM -0700, pcr...@pcreso.com wrote:
> Hi,
>
> I'm trying to determine all the possible intersections in sets of polygons.
>
> The approach I'm using is to load each set of polygon perimeters into a Postgis topology & extract all the faces as geometries. This results in some new polygons lying between the input ones. Which is the nature of this technique- that is not the problem.
>
> I then determine the number of original polygons that the new ones overlap, which should go from 0 up. My problem arises in that I'm getting new polygons, derived from the boundaries of the old ones, which are actually between old polygons, but have a non-zero (but very small - like 1e-15) area of overlap - but all they have in common is the boundary line which should have a zero area.
>
> So some of my new polygons supposedly overlap an old polygon, but in fact they do (or should) not.

How are you checking the overlap area ?
And how are you checking what polygons have "in common" ?

The thing is that building the topology implies noding the input
which in turn can change its shape. As a result the final edges are not
always fully covered by the original linework.

--strk;

[Brent Wood]

Hi,

Any advice appreciated!!

I'm undertaking a spatial analysis using Postgis (what else would I use!!!). The first part works well.

I take a large number (potentially millions) of lines defined by start & end points & buffer them to create polygons. (I'm working in lat/long EPSG:4326 but transforming to a custom equal area projection for the buffering operation).

I generate a grid of 5x5km cells (polygons) covering the region of interest.

I clip the line based polygons to the grid, so I can generate statistics for each cell describing the lines that intersect with it, various quantitative measures such as ST_Union() the clipped line polygons to generate a footprint in each cell to work out how much is/is not covered, or sum the ST_Area() of the clipped polygons grouped by cell to calculate an aggregate cover, which can be several times the actual cell area.

So far so good, it works well, the code is clear & transparent & provides a good result. At least as good as any commercial software can do. My test data subset is processed from scratch in about 30 minutes.

Now I want to ST_Union() all the cell based polygons into an overall single multipolygon representing the footprint. The code is simple. The performance, even with my subset,  is a problem.

I have thousands of cell based footprint multipolygons, each potentially with thousands of vertices to be ST_Union()ed. Runtime is weeks for an iteration. If I need separate total footprints for 20 different species annually for 5 years, that is 100 iterations. Memory & I/O use is minimal - it is totally cpu bound.

I am looking at trying to simplify the polygons to be unioned to reduce the number of vertices (& hence processing) involved, but to achieve any significant benefit I'm having to change the shape of the polygons to ST_Union() too much.



Does anyone have any suggestions as to how this could be made significantly faster?
If I had $$ to throw at developers to work on the codebase (presumably GEOS?) could performance be significantly improved?


Thanks,

   Brent Wood

[Bborie Park]

Your best bet is to consider splitting the workload among several postgresql connections.

darkblueb had a blog post about this...

http://blog.light42.com/wordpress/?p=23

_______________________________________________
postgis-devel mailing list
postgi...@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-devel

[map...@light42.com]

yes, true.. but with a thorough read you might notice that the gdal_retile.py experiment was largely ineffective, 

but if you click on the link at the top to the *next post*  

      Variable Buffers in PostGIS

you will find the one that really worked well.. in fact, we used that 2nd post in production for months, to great effect.

The trick on one machine was to split to work by some constant, and then make psycopg2 connections for each "bucket."

This worked very well..

Since then I have experimented only a tiny bit with SPARK from the Berkeley Amp Lab for a distributed work load on a Hadoop file system, but that world has no GEOS (yet) 

--

Brian M Hamlin

OSGeo California Chapter

blog.light42.com

[Brent Wood]

These tend to assume each operation is constrained to a single cell, hence parallelizable, hence undertaking an operation on multiple cells concurrently.

While I'm dealing with many cells - they are all being merged into a single multipolygon feature. I can't load two cells into the same multipolygon at the same time - two writes to the same record - so can't run concurrently.

What may be possible is to perhaps run multiple processes on subsets to create intermediate (larger) merged polygons which can then be merged themselves to create the final single feature.

This would probably allow better use of resources on a multi core system... at present I'm using 1.7% of memory on a 100% cpu process, so I'll look into this approach - 8 cores running concurrently giving giving close to 8x faster is useful.

Brent

---

From: "map...@light42.com" <map...@light42.com>
To: Brent Wood <pcr...@pcreso.com>; PostGIS Development Discussion <postgi...@lists.osgeo.org>; Bborie Park <dust...@gmail.com>
Cc: PostGIS Users Discussion <postgi...@lists.osgeo.org>
Sent: Thursday, October 17, 2013 3:16 PM
Subject: Re: [postgis-devel] ST_Union() performance problem (with possiblefunding)

[Rémi Cura]

Hey,

I'm not sure I understand perfectly your problem, 

for instance why exactly do you need an unioned footprint? 

Could you not stick to function where f(footprint) = sum(f(cells)) (this work for area, cover analysis etc.).

Meaning you could aggregate your result cell by cells.

Another think is what about your input data precision? I'm guessing it is low precision, so you may cut corners during processing (for instance, using buffer with quad_segs=2 to lower number of vertices in generated buffer, so generated footprint).

If not, and for very big data you may need to change your way of working.

The ideal candidat (lots of memory, less cpu) is image processing!

Fortunately you can switch to image world, where operations are way way faster, and easily made in parallel.

In image world you don't have to care about topology and exact computation, you can just compute by pixels.
(I won't even talk about GPU processing, but it may be used depending on your operations)

(If you don't require too much precision and can cope with aliasing effect)

Your data process would look like this :

using gdal directly (or several connection to postgis) you rasterize your polygons cells by cells. You do it the parallel way. You should have one image file per cell as output.

This images can be assembled into one big image depending on memory limitation.

If it doesn't fit in memory, you can use very fast , very efficient out of memory parallel processing software like Orfeo ToolBox, included in QGis sextante, depending on what you want to do.

What is cool with images is that you can do things you can't with pure geometry : 

I'm guessing you are working on data that are observations of animals. Instead of using a plain surface like "area where it is likely animals live", you could use fuzzy models (in the area, probability of animal living would be decreasing as the distance with the observation of animal),

then generate more "precise" and interesting results.

I made a lot's of guess so my answer may be completely useless .

I can precise things unclear.

Cheers,

Rémi-C

2013/10/17 Brent Wood <pcr...@pcreso.com>

[Brent Wood]

Hi,

I have a dataset with start & end points that we are represtnting as linestrings. Many records do not have endpoints - so we are takimg the median bearling for comparable startpoints & creating an arbitrary endpoint.

While it is easy to get a bearing for a known endpoint with ST_Azimuth(), is there a simple way in Postgis to get the point at a defined distance & bearing? I've perused the docs & nothing jumped out as an obvious solution.

Thanks,

  Brent Wood

---

From: Simon Greener <si...@spatialdbadvisor.com>
To: PostGIS Users Discussion <postgi...@lists.osgeo.org>
Sent: Thursday, April 11, 2013 6:57 PM
Subject: Re: [postgis-users] Finding closet intersect point along direction of line

[David Rowley]

On Sun, Dec 15, 2013 at 2:52 PM, Brent Wood <pcr...@pcreso.com> wrote:

Hi,

I have a dataset with start & end points that we are represtnting as linestrings. Many records do not have endpoints - so we are takimg the median bearling for comparable startpoints & creating an arbitrary endpoint.

While it is easy to get a bearing for a known endpoint with ST_Azimuth(), is there a simple way in Postgis to get the point at a defined distance & bearing? I've perused the docs & nothing jumped out as an obvious solution.

Sounds like you need ST_Project()

 

Regards

David Rowley

[Brent Wood]

Thanks David - I do indeed...

Cheers,
 

Brent

---

From: David Rowley <dgrow...@gmail.com>
To: Brent Wood <pcr...@pcreso.com>; PostGIS Users Discussion <postgi...@lists.osgeo.org>
Sent: Sunday, December 15, 2013 6:08 PM
Subject: Re: [postgis-users] How can I create a point at a known distance & bearing from a point in Postgis?