Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
1  
It looks like you're "doing it the hard way" in terms of the queries, but since you've already put so much work in to this design, maybe best not to change it. One option that keeps the current design might be to split your "shared features" - like roads, rivers, etc at the country border (so turn a single shared feature into several non-shared pieces). That might allow you to move forward with the design you've adopted. –  BradHards Dec 21 '12 at 11:00
    
thanks for the answer. However it is the type of query needed to do this. Ideally adding a country_name to each table would speed up the query, because I want to avoid having to split the feature class table up into separate country tables. –  Daily Oliver Dec 22 '12 at 14:08
    
Chris, so that you can comment on your own question, I have merged your two accounts. –  whuber Dec 22 '12 at 17:51
    
Please report the count of rows into your roads table: select count(*) from roads to know if it is too many roads the problem. Could you delete road outside country named 'Sudan'? Also how did you loaded the countries into database? Did you used the shp2pgsql to create sql file ? –  cavila Dec 28 '12 at 2:30
add comment

1 Answer

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;
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.