[GIS] Update a newly added column of a table using ST_Intersection as a filter

Question

I would be really grateful to anyone that could help me with the following issue.
I have ever used PostGIS much before – so could be considered a newbie…

My Situation:

I have TWO PostGIS tables:

1. roads
2. country_boundary

The roads table is so big (within the application 500MB) that I wish to display it by country_boundary.

I can do this successfully using the following query:

select ST_ASBINARY(roads.geog::geometry) from roads, country_boundary 
where country_boundary.adm0_name='Sudan' 
and ST_Intersects(roads.geog::geometry, country_boundary.geog::geometry)

HOWEVER it takes too long to display the data (several minutes) note: I have tried most optimization techniques, spatial indexing, (ST_Simplify/ST_SimplifyPreserveTopology) these either remove data or takes just as long.

I thought I had found a simple solution that worked – I add a new column 'adm0_name' (aka country name) to the roads table and populated it using the ST_Intersects query above.

update roads 
set adm0_name = 'Sudan'
where ST_ASBINARY(roads.geog::geometry) in
(
  select ST_ASBINARY(roads.geog::geometry) from roads, country_boundary where country_boundary.adm0_name='Sudan' and ST_Within(roads.geog::geometry, country_boundary.geog::geometry)
);

The resulting display query: works fast, and as expected.. >>problem solved<< ??

select ST_ASBINARY(roads.geog::geometry) from roads where roads.adm0_name='Sudan'

UNFORTUNATELY NO…

The ST_Intersects includes a small number of geometries (polyline(s) that are shared with neighbouring countries. This means polylines (roads) crossing the country boundary, cannot be give a unique country name. Note: I originally used ST_Within but this removed data that was shared between neighbouring countries (as expected with features like roads/rivers that do not respect country boundaries).

….. thanks for bearing with me…

I decided I needed to use ST_Intersection because this would avoid the boundary issue – it would 'clip' the road dataset to each country boundary thus avoiding the polyline boundary issue.

This works – I have an intersection query that displays the data OK

Select ST_ASBINARY (
ST_Intersection(g1.geom1, g1.geom2)) FROM (Select country_boundary.geog::geometry As geom1, roads.geog::geometry As geom2 from country_boundary, roads 
where country_boundary.adm0_name = 'Sudan') AS g1;

MY PROBLEM: I cannot find a way to update the roads.adm0_name column with the country name because ST_Intersection seems to returns a geometry, it does not behave like ST_Intersects, where I could update the column.

My attempt to do this is below —> 26 hours later I have reached a point of head scratching. Could anyone of you guys help me out?

update roads 
set adm0_name = 'Sudan'
where ST_ASBINARY(roads.geog::geometry) in
(
 Select ST_ASBINARY (
ST_Intersection(g1.geom1, g1.geom2)) FROM (Select country_boundary.geog::geometry As geom1, roads.geog::geometry As geom2 from country_boundary, roads where country_boundary.adm0_name = 'Sudan') AS g1;
);

However the above statement fails —> I need to update the roads.adm0_name with country_boundary.adm0_name where (country_boundary "ST_Intersection" roads).

Thanks again —

Answer 1

Thanks for the help.

In the end I made a new table and copied over the columns I needed -- solution below Note: worthy mention to Paul Ramsey who recommended the method below to quicken up the ST_Intersection the query went from 6 hours to 4 minutes. Worth looking up...

BEGIN;

--insert into contours

insert into contours_country(gid,fnode_,tnode_,lpoly_,rpoly_,length,contourl_,contourl_i,f_code,f_code_des,hqc,hqc_descri,zv2,zv2_descri,tile_id,edg_id,adm0_name,geog,geom)
select 
a.gid As gid,
a.fnode_ As fnode_,
a.tnode_ As tnode_,
a.lpoly_  As lpoly_,
a.rpoly_ As rpoly_,
a.length As length,
a.contourl_ As contourl,
a.contourl_i AS contourl_i,
a.f_code As f_code,
a.f_code_des As f_code_des,
a.hqc As hqc,
a.hqc_descri As hqc_descri,
a.zv2 As zv2,
a.zv2_descri As zv2_descri,
a.tile_id As tile_id,
a.edg_id As edg_id,
b.adm0_name As adm0_name,
a.geog As geog,

CASE
WHEN ST_WITHIN(a.geog::geometry,b.geog::geometry)
THEN a.geog::geometry
ELSE ST_Intersection(a.geog::geometry, b.geog::geometry)
END AS geom

From contours a
JOIN country_boundary b

ON ST_Intersects(a.geog::geometry, b.geog::geometry)

Where b.adm0_name ='Niger';

--insert in roads

insert into roads_country(gid,fnode_,tnode_,lpoly_,rpoly_,length,road3a_,road3a_id,fnode___7,tnode___8,lpoly___9,rpoly___10,road3_,road3_id,rdline_,rdline_id,rdlntype,rdlnstat,adm0_name,geog,geom)
select 
a.gid As gid,
a.fnode_ As fnode_,
a.tnode_ As tnode_,
a.lpoly_ As lpoly_,
a.rpoly_ As rpoly_,
a.length As length,
road3a_ As road3a_,
road3a_id As road3a_id,
fnode___7 As fnode___7,
tnode___8 As tnode___8,
lpoly___9 As lpoly___9,
rpoly___10 As rpoly___10,
road3_ As road3_,
road3_id As road3_id,
rdline_ As rdline_,
rdline_id As rdline_id,
rdlntype As rdlntype,
rdlnstat As rdlnstat,
b.adm0_name As adm0_name,
a.geog As geog,

CASE
WHEN ST_WITHIN(a.geog::geometry,b.geog::geometry)
THEN a.geog::geometry
ELSE ST_Intersection(a.geog::geometry, b.geog::geometry)
END AS geom

From roads a
JOIN country_boundary b

ON ST_Intersects(a.geog::geometry, b.geog::geometry)

Where b.adm0_name ='Niger';



--insert into wadis

insert into wadis_country(gid,fnode_,tnode_,lpoly_,rpoly_,length,dnnet_,dnnet_id,dnlntype,dnlnstat,adm0_name,geog,geom)
select 
a.gid As gid,
a.fnode_ As fnode_,
a.tnode_ As tnode_,
a.lpoly_ As lpoly_,
a.rpoly_ As rpoly_,
a.length As length,
a.dnnet_ As dnnet_,
a.dnnet_id As dnnet_id,
a.dnlntype As dnlntype,
a.dnlnstat As dnlnstat,
b.adm0_name As adm0_name,
a.geog As geog,

CASE
WHEN ST_WITHIN(a.geog::geometry,b.geog::geometry)
THEN a.geog::geometry
ELSE ST_Intersection(a.geog::geometry, b.geog::geometry)
END AS geom

From wadis a
JOIN country_boundary b

ON ST_Intersects(a.geog::geometry, b.geog::geometry)

Where b.adm0_name ='Niger';

COMMIT;