[postgis-users] Discovering where roads cross rivers

[Stephen Woodbridge]

Dave Potts wrote:
> I am trying to discover all the data points when a river system crosses 
> a road system, my data is expressed as two different shapefiles,  I had 
> assumed that it was just a case of calling  the Intersection method
> 
> eg. select m.name from transport m , rivers r where 
> Intersection(r.the_geom,m.the_geom) = 't';
> 
> But I get a 'parse error- invaid geometry when I try this, is this the 
> correct way of doing this?

Start by:

select count(*) from transport where not is_valid(the_geom);
select count(*) from transport where not is_valid(the_geom) and not 
is_valid(buffer(the_geom, 0.0));

If the later clears up all the issues then do:

update transport set the_geom = buffer(the_geom, 0.0);

and do the same for rivers.

If you still have errors you will need to load those geometries into 
other software to fix them. of you can ignore them by filtering them out 
of the intersection query.

Also intersection returns the geometry of the intersection. if you just 
want the boolean it is intersects() and you should add a gist index to 
the_geom column on both tables and

select m.name from transport m , rivers r where
   m.the_geom && r.the_geom and                       -- use index
   is_valid(m.the_geom) and is_valid(r.the_geom) and  -- add filter bad
   intersects(r.the_geom,m.the_geom);                 -- boolean test

-Steve W