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 have a PostGis database created with osm2pgsql. I am looking for a way to add the country name in which a point is located to the point row. It must work also for countries with enclaves, such as France (Llívia). In this case, osm2pgsql creates several polygons from the border relation: a polygon with the outer border and and the enclave, and for each of the both countries a polygon with the enclave only. A simple ST_Within or ST_Contains does not work in this case. I found I way to do the job:

update planet_osm_point p set country=
       (SELECT distinct poly.name 
         from planet_osm_point point inner join planet_osm_polygon poly 
         on (st_within(point.way, poly.way) 
            and (st_nrings(poly.way)=1 or not st_within(point.way, st_makepolygon(ST_InteriorRingN(poly.way,1))))
            and not exists (select * from planet_osm_polygon other 
                              where other.admin_level='2' 
                                    and other.name=poly.name 
                                    and st_contains(st_makepolygon(st_exteriorring(other.way)), poly.way)
                                    and other.way_area>poly.way_area))
            where  poly.boundary='administrative' and poly.admin_level='2' and point.osm_id=p.osm_id) 
      where p.name is not null;

This seems a bit complex to me and takes a lot of processing time. Is there any simpler and faster way to do this?

share|improve this question
    
I found the mistake that I made. I merged region extracts and country border extracts with osmosis and then imported with osm2pgsql to the database. This gave me an extra boundary polygon for the enclave and the surrounding country. It was this polygon which caused problems, as st_within returned true when applied to a way from the enclave. I now import the country border in a second osm2pgsql run and everything works fine. –  user17607 May 1 '13 at 7:02
add comment

1 Answer

If the enclaves from the polygon country layer are true "donut holes" then ST_Intersects() should work as you expect. i.e. (you probably already have a spatial index on the countries?)

update planet_osm_point p set country=
       (SELECT distinct poly.name 
         from planet_osm_point point inner join planet_osm_polygon poly 
         on (st_intersects(point.way, poly.way) 
         where  poly.boundary='administrative' and poly.admin_level='2' and point.osm_id=p.osm_id) 
      where p.name is not null;

If the country polygons do not have true holes for the enclaves, then your (clever) workaround is probably what's needed...

share|improve this answer
    
osm2psql creates a polygon with an outer ring and an inner ring for the enclave. So there is a true hole and st_intersects works well for a way in the enclave and this polygon. –  user17607 May 1 '13 at 6:51
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.