0

I have a Postgres database with a postgis extention installed and filles with open street map data.

With the following SQL statement :

SELECT                                                                
    l.osm_id,
    sum(
        st_area(st_intersection(ST_Buffer(l.way, 30), p.way))
        /
        st_area(ST_Buffer(l.way, 30))
    ) as green_fraction
FROM planet_osm_line AS l
INNER JOIN planet_osm_polygon AS p ON ST_Intersects(l.way, ST_Buffer(p.way,30))
WHERE p.natural in ('water') or p.landuse in ('forest') GROUP BY l.osm_id;

I calculate a "green" score.

My goal is to create a "green" score for each osm_id.

Which means; how much of a road is near a water, forrest or something similar.

For example a road that is enclosed by a park would have a score of 1.

A road that only runs by a river for a short period of time would have a score of for example 0.4

OR so is my expectation.

But by inspection the result of this calculation I get sometimes Values of

212.11701212511463 for a road with the OSM ID -647522

and 82 for a road with osm ID -6497265

I do get values between 0 and 1 too but I don't understand why I do also get such huge values.

What am I missing ?

I was expecting values between 1 and 0.

Andreas
  • 397
  • 4
  • 18
  • 37
  • 3
    Duplicate osm_ids that you are summing. Try without summing and grouping, to see the duplicate results for the IDs you mention. – pkExec Jan 18 '23 at 10:34
  • the calculation takes around 4 days so I can only see if thats correct in 4 days :) But it seems to me your suggestion is plausible. Thank you. – Andreas Jan 18 '23 at 12:22
  • 1
    if you do a quick test using for example "WHERE l.osm_id=-6497265 AND (p.natural........ OR .......)", you should get an instant response, provided that you have an index on osm_id, and a spatial index on the "way" column. – pkExec Jan 18 '23 at 12:53
  • Another source of issue is that multiple polygons could overlap each others. You would have to `st_union` their geometries first. – JGH Jan 18 '23 at 14:05
  • On a side note, you can join using `st_dwithin()` instead of buffer+intersection. It will make use of the index and will be much faster. – JGH Jan 18 '23 at 14:06
  • could you please provide an example on how would I use st_dwithin() and st_union()? – Andreas Jan 18 '23 at 15:03
  • This was mine approach: SELECT l.osm_id, st_area(ST_DWithin(l.way, p.way,30)) / st_area(ST_Buffer(l.way, 30)) as green_fraction FROM planet_osm_line AS l INNER JOIN planet_osm_polygon AS p ON ST_DWithin(l.way,p.way,30) WHERE l.osm_id=-6497265 AND p.natural in ('water') or p.landuse in ('forest'); and I get No function matches the given name and argument types. You might need to add explicit type casts. – Andreas Jan 18 '23 at 15:09

1 Answers1

1

Using a custom unique ID that you must populate, the query can also union eventually overlapping polygons:

SELECT                                                                
    l.uid,
    st_area(
        ST_UNION(
            st_intersection(ST_Buffer(l.way, 30), p.way))
    ) / st_area(ST_Buffer(l.way, 30)) as green_fraction
FROM planet_osm_line AS l
INNER JOIN planet_osm_polygon AS p 
    ON st_dwithin(l.way, p.way,30)
WHERE p.natural in ('water') or p.landuse in ('forest') 
GROUP BY l.uid;
JGH
  • 15,928
  • 4
  • 31
  • 48
  • just a follow up:this query worked beautifully. The previous query took around 4 days to complete. The suggested one by JGH took around 6 hours on my machine. – Andreas Jan 19 '23 at 07:41