How to find the country in which a point is located?

Question

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?

Answer 1

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...

Answer 2

I solved this by adding this subselect to my query:

select name where admin_level = 2 and (select geometry from place where
  osm_id=p.osm_id) && geometry and ST_Within((select geometry from place where
  osm_id=p.osm_id), geometry) limit 1;

This is a query against the gazetteer database but you could translate this to use against the standard osm schema.

I could not get ST_Within to use the spatial index for some reason so I added the intersects clause (&&) to help out.