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 am in the process of building a data warehouse in PostGIS, that will store the length of the road network of a country (Lines) per Administrative Area (i.e. Province or State, referred to as AA hereafter and is a Polygon). This is my first time of building a large scale Relational data warehouse so I am not really familiar with this kind of operations.

The "id" field of the spatial join resulting table is used to relate to the unique "id" field of the AA polygon layer. The length of the road segments is already calculated and is derived from the "len" attribute of the Road Network table.
So far the operation yields the desired results by performed the following query:

SELECT  
  m.id,  
  m.order08,  
  sum(r.len) as length,  
FROM
  public.and_lines_2012_03 AS r, 
  public.and_a8_2010_12 AS m
WHERE
  ST_Intersects(m.geom,r.geom) 
GROUP BY 
  m.id, m.order08
ORDER BY 
  id;

And the resulting table looks something like this:

gid |   id   | length  
 1  | 000002 | 74118
 2  | 000012 | 128131
 3  | 000013 | 174296
 4  | 000014 | 82240
 ...|  ...   |  ...

The Road Network is divided in 3 different line types that are attributed in the Road Network table by the name "line_type" and can have values 0,1,2 My intention is to create a table that will have a column for the length of each line type for each Administrative Area. That is the table would ideally look something like this:

gid |   id   | length_type_00 | length_type_01 | length_type_02
 1  | 000002 |     30028      |     11564      |     32526
 2  | 000012 |     64270      |     15602      |     48259
 3  | 000013 |     104492     |     51981      |     17823
 4  | 000014 |     38153      |     9524       |     34563
 ...|   ...  |      ...       |      ...       |      ...

I have experimented with different approaches in SQL statements, but so far have not managed to grasp the approach of this kind of queries. The most relative resource I have found on the Internet is this one from OpenGeo:

Any help would be highly appreciated

Kind regards
Sergios

Update 2012/05/23

@Nicklas

I have managed to produce what seems to be very accurate results with the way you have suggested. The query used is the following:

SELECT  
  b.id AS id,
  b.order08 AS order08,
  round(sum(b.line_type_00)) AS line_type_00,
  round(sum(b.line_type_01)) AS line_type_01,
  round(sum(b.line_type_02)) AS line_type_02
FROM
(SELECT  
  m.id,
  m.order08,
  (line_type=0)::int * sum(ST_Length(ST_Intersection(r.geom::geography, m.geom::geography))) line_type_00,  
  (line_type=1)::int * sum(ST_Length(ST_Intersection(r.geom::geography, m.geom::geography))) line_type_01,  
  (line_type=2)::int * sum(ST_Length(ST_Intersection(r.geom::geography, m.geom::geography))) line_type_02
FROM
  public.and_lines_2012_03 AS r, 
  public.and_a8_2010_12 AS m
WHERE
  ST_Intersects(m.geom,r.geom)
GROUP BY 
  m.id, m.order08,line_type
ORDER BY
  m.id
) 
AS
  b
GROUP BY
  b.id,b.order08
;  

The query is quite slow though (4846 tuples take around 41 seconds), since from what I understand it re-projects on the fly to produce the the results in meters. Have I understood/am I doing something wrong?

Thank you in advance
Sergios

share|improve this question
    
In what srid is your data? Your cast to geography doesn't make sense to me. –  Nicklas Avén May 23 '12 at 15:47
    
@NicklasAvén The SRID is 4326 (WGS_1984) imported from GCS_1984 shapefile –  Sergios Kolios May 23 '12 at 15:55
    
Ok, sorry. I didn't know that ST_Intersection supports geography type, but it does. And as you say it is projected on the fly. If it is possible with respect to the extent of your data, you would gain a lot of projecting the data to a suitable projection before doing the calculations. –  Nicklas Avén May 23 '12 at 20:00
    
@NicklasAvén Thank you for all the help. I will probably do as you suggest, that is map projected coordinate systems per country in a table of the database and do the projection upon import, I think this will be computationally profitable. The number of columns for my example is fixed, so your reply is actually the reply to my question, but I will leave the question for a few days in case somebody gives a reply that includes pivot-tables. Thanks again –  Sergios Kolios May 23 '12 at 20:19
    
Here is the doc for pivot functionality. postgresql.org/docs/9.1/static/tablefunc.html. But I think it is a less transparent way of doing it. But maybe it is faster. I don't know. –  Nicklas Avén May 23 '12 at 20:27
add comment

1 Answer 1

up vote 4 down vote accepted

What you want is a pivot table. You can create that with the pivot functionality in PostgreSQL, but I think it is easier to create it like this in your case:

SELECT  
  m.id,  
  (length_type='00')::int * sum(ST_Length(r.geom)) length_type_00,  
  (length_type='01')::int * sum(ST_Length(r.geom)) length_type_01,  
  (length_type='02')::int * sum(ST_Length(r.geom)) length_type_02
FROM
  public.and_lines_2012_03 AS r, 
  public.and_a8_2010_12 AS m
WHERE
  ST_Intersects(m.geom,r.geom) 
GROUP BY 
  m.id, m.order08
;

But be aware here that you will get to many km of roads. The roads intersecting two polygons will be represented in both with their full length added to both polygons. Depending on your data you have two options to solve this.

I your roads stops at the border to the polygon, it should, in it's full length be added to one of the polygons. Then you can use ST_Line_Interpolate_Point(geom, 0.5); to get the middle of the road. Then use that point in your ST_Intersection. Then the road will only be joined to the polygon with the roads center point.

The second option is if your roads crosses the border. Then, to get the right result you will have to use ST_Intersection to get the portion of the road that is in that particular polygon. That will be slower if you have a lot of roads and or the roads or polygons have many vertex-points.

The query could then look:

SELECT  
  m.id,  
  (length_type='00')::int * sum(ST_Length(ST_Intersection(r.geom, m.geom))) length_type_00,  
  (length_type='01')::int * sum(ST_Length(ST_Intersection(r.geom, m.geom))) length_type_01,  
  (length_type='02')::int * sum(ST_Length(ST_Intersection(r.geom, m.geom))) length_type_02 
FROM
  public.and_lines_2012_03 AS r, 
  public.and_a8_2010_12 AS m
WHERE
  ST_Intersects(m.geom,r.geom) 
GROUP BY 
  m.id, m.order08
;

HTH

Nicklas

share|improve this answer
    
Hi Nicklas Thank you very much for the answer, indeed I am looking for pivot-table functionality. –  Sergios Kolios May 23 '12 at 10:45
    
I have realised that I get "too many roads" when comparing the sum of stats for all AA's with the total stats of the Country, significant deviations occur in relation to the expected numbers. I will attempt to implement it according to your suggestions and post back the results. Thank you for the effort Sergios –  Sergios Kolios May 23 '12 at 11:37
    
Yes, but remember that you then have to recalculate the lengthes as in my example. The pre calculated will not help you if you use ST_Intersection to get the share of the road in each polygon. But at the other side, ST_Length is quite, very, fast :-) –  Nicklas Avén May 23 '12 at 12:33
    
I think I have done what you have suggested and came up with pretty accurate results. You can find above the query used. The only issue now is speed, since I will have to work with extremely large datasets and on the fly projection seems to take a big toll on the performance. –  Sergios Kolios May 23 '12 at 14:59
    
@Sergios Kolios Ops, I see that, of course it is necessary to divide the query like you have done in two levels. I will update some day –  Nicklas Avén Jun 1 '12 at 11:54
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.