Calculating length of just portion of linestring that intersect with polygon using PostGIS?

Question

I want to calculate the portion of the linestring that lies within the limits of a polygon. Better with an image

What is the length of the portion of linestring within Province A? and within Province B? I'm looking for a query that returns something like:

province_name | length   |
--------------------------
Province A    | 30.5342  |
Province B    | 19.4321  |

My approach is a PL/pgSQL script that loops over all the points in the linestring and set a mark in the first point that is in province B (point N). Then, calculate the length from point 0 to point N (distance within province A) and the length from point N+1 to last point (distance within province B). But sounds rusty, and difficult to extend to an indefinite number of provinces.

Is there any better approach? Suggestions welcome. I'm using PostGIS 2.0, by the way.

---

@HeyOverThere suggests a query that almost works, with a little modification:

SELECT a.name, sum(ST_LENGTH(ST_Intersection(a.geom, b.geom))/1000)
FROM polygons a, lines b
WHERE ST_Intersects(a.geom, b.geom) group by a.name order by a.name

That query returns 2 columns, with the province name and a length. But the length is exactly twice the value it should be. Does it make sense?

More info: the provinces are really multipolygons (but all share a the province name as common field), and the line has been built from gps tracks (points) and st_makeline.

---

As @GetSpatial suggested, I've tried using subqueries to sum and group the results of the first query. Same result. The length of each intersection is correct. So, maybe the intersection is two times as long as I expect (?).

Here, the 3 things I've tried. All of them get a distance 2x the expected one (link to real data at the end):

First try, with one query:

SELECT a.province, sum(ST_LENGTH(ST_Intersection(a.geom, b.the_geom))/1000) as length_per_province
FROM provinces a, line_tracks b
WHERE ST_Intersects(a.geom, b.the_geom) group by a.province order by a.province;

Second try, with CTE:

with line_segments as (
SELECT a.province, ST_LENGTH(ST_Intersection(a.geom, b.the_geom))/1000 as length
FROM provinces a, line_tracks b
WHERE ST_Intersects(a.geom, b.the_geom) order by a.province
) select province, sum(length) as length_per_province from line_segments group by province;

Third one, with subquery, as suggested:

select province, sum(length) as length_per_province from(
SELECT a.province, ST_LENGTH(ST_Intersection(a.geom, b.the_geom))/1000 as length
FROM provinces a, line_tracks b
WHERE ST_Intersects(a.geom, b.the_geom)) as foo group by province;

I get this, in all cases:

"province";"total_length"
"p1";62.8687677452213
"p2";40.254207910339

The total length of the line:

select sum(st_length(the_geom))/1000 as total_length from line_tracks

The result

"total_length"
51.5614878277802

And the 2 tables here (SET client_encoding = 'LATIN1'): https://dl.dropboxusercontent.com/u/6599273/gis_data/data.sql.tar.bz2

Answer 1

Here's a query that should get you what you want:

SELECT a.name, b.name, ST_LENGTH(ST_Intersection(a.geom, b.geom))
FROM polygons a, lines b
WHERE ST_Intersects(a.geom, b.geom);