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

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

• I think in the process of add the group by set to your query, you are causing a double return of the values. Try taking the query that @HeyOverThere gave you, and wrapping that with another query that contains the grouping function. This way, you are doing the grouping and sum on the results of the first query, instead of as part of the first query. – Get Spatial Nov 20 '13 at 2:36
• Mmm... same problems. I've added my real data and the queries I've executed. So, if you want to try... – Jorge Arévalo Nov 20 '13 at 8:55
• @jorgeas80 have you figured this problem out? i am having a very similar problem – ziggy Nov 18 '16 at 16:53
• Not really. I had to deal with a different problem and didn't continue with this one (It was just a test). But if it's interesting for someone, I'm gonna give it a new try. Are you stuck in the same point? – Jorge Arévalo Nov 20 '16 at 11:29