Geographic city boundaries

Hi all!

I’m getting the city boundaries with this query:

        (SELECT * FROM
          (SELECT
              planet_osm_polygon.way,
              planet_osm_polygon.admin_level,
              planet_osm_point.place,
              CASE
                WHEN (planet_osm_point.tags->'population' ~ '^[0-9]{1,8}$') THEN (planet_osm_point.tags->'population')::INTEGER ELSE 0
              END as population,
              (ST_Area(ST_Transform(planet_osm_polygon.way, 4326)::geography) / 1000000) AS km2
            FROM planet_osm_polygon
            JOIN (
                   SELECT name, MAX(admin_level) AS al
                   FROM planet_osm_polygon
                   WHERE boundary = 'administrative' AND admin_level IN ('4', '6', '8') AND osm_id < 0 GROUP BY name
                 ) size USING(name)
            JOIN planet_osm_point USING (name)
            WHERE planet_osm_polygon.boundary = 'administrative' AND
                  planet_osm_polygon.admin_level = size.al AND
                  (
                    (
                      planet_osm_polygon.admin_level IN ('6', '8') AND
                      planet_osm_point.place IN ('city', 'town')
                    ) OR
                    (
                      planet_osm_polygon.admin_level = '4' AND
                      planet_osm_point.place = 'city'
                    )
                  ) AND
                  planet_osm_polygon.osm_id < 0 AND
                  planet_osm_polygon.way_area > 1*!pixel_width!::real*!pixel_height!::real
          ) AS boundaries
          WHERE (admin_level = '4' AND km2 >= 100) OR
                (
                  population > 20000 AND
                  (
                    (admin_level = '6' AND km2 >= 50) OR
                    (admin_level = '8' AND km2 >= 25)
                  )
                ) OR
                (
                  population > 10000 AND admin_level = '8' AND km2 >= 100
                )
        ) AS city_boundaries

So I can get the political boundaries (administrative) of cities the have some requisites.
It works, but for my maps (using by flying) the geographical boundaries are better. I can see many little town in my map with a big area, but on the area is just “pampa”.
Is it possible to restrict the area to the part of the city where the are more than X buildings?

Thanks a lot
Luca

I did some experimental work to derive urban/built-up areas from OSM. You can see my OSM-based data here. This will be out-of-date because I did it in 2016. The data is available as geojosn tiles of 10 square degrees.

Thank you for your answer!
I’d like to try with Saxony, since I’m here and I have many maps to compare the work…
Could you explain me how you did this work?

Thanks
Luca

Maybe a clarification…
With “geographic city boundaries” I mean: the part of the city where I can see buildings (not just an isolated house, but a collection of buildings).

Thanks for your help
Luca

You can use ST_ClusterDBSCAN() and ST_ConcaveHull().

Thank you very much for your suggestion, but… on which data? :stuck_out_tongue:
Could you suggest me a query to find the needed data?

Thanks
Luca

I would think on buildings, but the approach using residential roads or short length roads would work in similar ways. The advantage of using roads is this works even if no buildings are mapped. Some positive & negative buffering might be needed to fill holes.

OK, but could you explain me how to do that?
I really don’t have any idea how to create a polygon with these data…

Thanks
Luca

So, maybe I found a way, but I need more help for the implementation…
I can use the landuse = ‘residential’, with a query like that:

(SELECT
               way
           FROM planet_osm_polygon
           WHERE landuse IN ('residential', 'retail', 'retail;residential')

So I can select the areas with building, that is what I want…
My problem is, that I want to avoid isolated buildings and create a boundary of the “main area”, to avoid having a map like this (residential are the red spots):

Can someone suggest me a way?

Thanks
Luca

The problem of landuse areas is a) they are not always mapped; b) you need a few others commercial, industrial; c) you also need a few insitutions schools, universities, hospitals; d) some of these landuses may be ex-urban.

When I get time today I’ll show you a SQL query for clustering buildings.

Hi

Well, I can add the other types, this will not the problem… :wink:
My problem is to have a boundary that I can measure, so that I can decide to only display areas greater than X. If I can get that, ex-urban and little town are filtered by default…

Thanks
Luca

Hi again,

Maybe I’m on the right way, but I don’t know how to prosecute my idea…
I see the function ST_ClusterDBSCAN, so (if I understand correctly!) I can have an “ID” of contiguos elements.
So I wrote this query:

SELECT way, ST_ClusterDBSCAN(way, eps := 50, minpoints := 2) OVER () AS cid FROM planet_osm_polygon WHERE landuse = 'residential'

Nice… And now? How can I create a geometry containing all elements grouped by “cid”?

Thanks a lot for your help
Luca

Sorry took me a bit longer, here’s an example based on a local authority (but you can change the st_within to a bbox:

select cid, st_buffer(st_buffer(st_concavehull(st_multi(st_union(b)),0.9),100),-100) as geom, count(*) as build_cnt from (
select row_number() over() as id, osm_id, st_clusterdbscan(way,100,5) over() as cid, st_centroid(way) as b
from planet_osm_polygon
where building is not null
and st_within(way,(select way from planet_osm_polygon where name =‘Nottinghamshire’ and admin_level=‘6’))) a
where cid is not null
group by cid

I use a minimum of 5 buildings within 100 m (in pseudomercator) for clustering. It’s important to filter out stuff which fails to be clustered (cid is not null). Union the points or ways which ought to create a multipolygon or multipoing, but force it to multi to be sure, then buffer by 100 m & re-buffer by the negative amount (which should fill small holes).

You may need to change the st_clusterdbscan parameters. 5 buildings seems to find a lot of farms, so probably should be higher. 100-150 m generally works well on UK data.

You are an hero!!!
I tried your query (just leaving the condition for Nottinghamshire) and it seems good… I just have to filter out all little towns, but I can do that moving your query in a subquery and adding WHERE (ST_Area(ST_Transform(way, 4326)::geography) / 1000000) > 40

Only disadvantage: your query took about half an hour to render a little part of Saxony, but it can be acceptable since I already know I need months to render Europe…
Maybe I can fire the query and save the result in a table to be faster in all rendering levels? Can I just run

CREATE TABLE city_boundaries AS ...

with your query and then in Mapnik

SELECT way FROM city_boundary WHERE (ST_Area(ST_Transform(way, 4326)::geography) / 1000000) > 40

?

I’ll post the result here as soon I got the “filter for little towns” working.

Thanks
Luca

The slowest bit of the query is probably the (select … where name=‘N…’) which I suspect does a scan on planet_osm_polygon. You can try replacing it with “where osm_id = ####” as this should use the primary index whereas AFAIK name is not indexed on planet_osm_polygon. Similarly an index on the building column may improve performance. I was more concerned to check why I got one big polygon as well as lots of little ones which is why the is null condition is needed. Agree with st_area, but also try upping the 5 to 10 or 20.

Also look at a materialised view as a way of storing the data in such a way that you can refresh it.

So, it does not seems that “where name = ‘xxx’” is the problem, but this is not important.
I’m trying to create a table with the data to use later in my queries, but I’d like to have the city name in the table. Of course I cannot use “name”, since in this case I’ll get the name of the building.
Could you suggest me a way to get the city name of a part of it (aka: I know the building and I want to know in which city is this building)?

Thanks
Luca

Not so easy (the urban areas will cover multiple named places).

You may want to rank cities by population & use the most populous in an area, or possibly use something like Voronoi triangulation to cut them up. There’s an example on my old blog post: http://sk53-osm.blogspot.com/2015/10/urban-areas-meditation-on-why-simple.html.

I’m trying to work another solution to my problem (reduce the “little spots” outside the cities) but I have a problem and I can’t find my error…

I created a table of the administrative city’s boundaries:

CREATE TABLE city_admin_boundaries AS
SELECT
    planet_osm_polygon.way,
    planet_osm_polygon.admin_level,
    planet_osm_point.name,
    planet_osm_point.place,
    CASE
     WHEN (planet_osm_point.tags->'population' ~ '^[0-9]{1,8}$') THEN (planet_osm_point.tags->'population')::INTEGER ELSE 0
    END as population,
    (ST_Area(ST_Transform(planet_osm_polygon.way, 4326)::geography) / 1000000) AS km2
  FROM planet_osm_polygon
  JOIN (
    SELECT name, MAX(admin_level) AS al
    FROM planet_osm_polygon
    WHERE boundary = 'administrative' AND admin_level IN ('4', '6', '8') AND osm_id < 0 GROUP BY name
  ) size USING(name)
  JOIN planet_osm_point USING (name)
WHERE planet_osm_polygon.boundary = 'administrative' AND
planet_osm_polygon.admin_level = size.al AND
(
  (
    planet_osm_polygon.admin_level IN ('6', '8') AND
    planet_osm_point.place IN ('city', 'town')
  ) OR
  (
    planet_osm_polygon.admin_level = '4' AND
    planet_osm_point.place = 'city'
  )
) AND
planet_osm_polygon.osm_id < 0;

and I have a table of the “buildings”:

CREATE TABLE city_boundaries AS
  SELECT NEXTVAL('city_boundaries_seq') AS id,
        ST_CollectionExtract(unnest(ST_ClusterWithin(way, 200)), 3)::geometry(MultiPolygon, 3857) as way
    FROM planet_osm_polygon
   WHERE landuse IN ('residential', 'retail', 'retail;residential', 'commercial', 'school', 'university', 'industrial',
                     'asphalt', 'cemetery', 'civic', 'civic_admin', 'concrete_surface', 'construction', 'education',
                     'educational', 'institutional', 'village');
ALTER TABLE city_boundaries OWNER TO _renderd;

-- This deletes all polygons which after removing 50 meters from all sides have area less than 50
DELETE FROM city_boundaries WHERE ST_Area(ST_Buffer(way, -50)) < 50;

-- Aggregates all clustered polygons by adding 120m and then removing 250m from all sides
UPDATE city_boundaries SET way = ST_Makevalid(ST_Multi(ST_Buffer(ST_Buffer(way, 250, 'join=miter'), -250, 'join=miter')));

With the buildings’ table I can draw what I need (I must surely adjust the polygons, but the way seems correct), but I have many buildings outsides the cities, so I try to get only the ones insides cities greater than X and with a population at least Y:

(SELECT
               city_boundaries.way
           FROM city_boundaries, city_admin_boundaries
           WHERE ST_Within(city_boundaries.way, city_admin_boundaries.way) AND
                 (city_admin_boundaries.population >= 18000 OR city_admin_boundaries.km2 >= 50)
          ) AS test

But so I have anyway many buildings outsides the cities and just few buildings inside the cities…

Can someone help me to find my error?

Thanks
Luca

Hi again,

So, I think I got a right way using a couple of help tables.
Now I can draw the buildings so (buildings boundaries are in red. Yellow are the administrative boundaries of the towns):
Since I can’t get the “holes” filled up in the database, I’m asking if it should be possible using CartoCSS…
Currently I have this code to select data:


  - id: TEST
    geometry: linestring
    <<: *extents
    Datasource:
      <<: *osm2pgsql
      table: |-
          (SELECT
               way
           FROM city_boundaries
           WHERE (ST_Area(ST_Transform(way, 4326)::geography) / 1000000) > 3.5
          ) AS test
    properties:
      cache-features: true
      minzoom: 7

and this one to draw them on the map:


#TEST {
  background/line-join: bevel;
  background/line-width: @city-boundaries-line-width;
  line-join: bevel;
  line-width: 2.5;
  line-color: #0000ff;
  background/line-color: #ff0000;
  polygon-fill: #ff0000;
  opacity: @city-boundaries-opacity;
}

Is there a settings in the CSS that I can use to fill up the holes?

Thanks
Luca

Hi!

I finally got it working!
Now it looks like that:

Pretty like the ICAO maps:

What I did:

  1. I created a help table with the administrative city boundaries:

CREATE TABLE city_admin_boundaries AS
SELECT
    NEXTVAL('city_admin_boundaries_seq') AS id,
    planet_osm_polygon.way,
    planet_osm_polygon.admin_level,
    planet_osm_point.name,
    planet_osm_point.place,
    CASE
     WHEN (planet_osm_point.tags->'population' ~ '^[0-9]{1,8}$') THEN (planet_osm_point.tags->'population')::INTEGER ELSE 0
    END as population,
    (ST_Area(ST_Transform(planet_osm_polygon.way, 4326)::geography) / 1000000) AS km2
  FROM planet_osm_polygon
  JOIN (
    SELECT name, MAX(admin_level) AS al
    FROM planet_osm_polygon
    WHERE boundary = 'administrative' AND admin_level IN ('4', '6', '8') AND osm_id < 0 GROUP BY name
  ) size USING(name)
  JOIN planet_osm_point USING (name)
WHERE planet_osm_polygon.boundary = 'administrative' AND
planet_osm_polygon.admin_level = size.al AND
(
  (
    planet_osm_polygon.admin_level IN ('6', '8') AND
    planet_osm_point.place IN ('city', 'town')
  ) OR
  (
    planet_osm_polygon.admin_level = '4' AND
    planet_osm_point.place = 'city'
  )
) AND
planet_osm_polygon.osm_id < 0;

  1. I create a second help table, with all the “main buildings” in the cities/towns starting from a given size:

CREATE TABLE city_boundaries AS
   SELECT NEXTVAL('city_boundaries_seq') AS id,
          ST_CollectionExtract(unnest(ST_ClusterWithin(planet_osm_polygon.way, 200)), 3)::geometry(MultiPolygon, 3857) as way
   FROM planet_osm_polygon, city_admin_boundaries
   WHERE landuse IN ('residential', 'retail', 'retail;residential', 'commercial', 'school', 'university', 'industrial',
                     'asphalt', 'cemetery', 'civic', 'civic_admin', 'concrete_surface', 'construction', 'education',
                     'educational', 'institutional', 'village') AND
         ST_Within(planet_osm_polygon.way, city_admin_boundaries.way) AND
            (
              (city_admin_boundaries.admin_level = '4' AND km2 >= 100) OR
              (
                population > 20000 AND
                (
                  (city_admin_boundaries.admin_level = '6' AND km2 >= 50) OR
                  (city_admin_boundaries.admin_level = '8' AND km2 >= 25)
                )
              ) OR
              (
                population > 10000 AND city_admin_boundaries.admin_level = '8' AND km2 >= 100
              )
            );
DELETE FROM city_boundaries WHERE ST_Area(ST_Buffer(way, -50)) < 50;
UPDATE city_boundaries SET way = ST_Makevalid(ST_Multi(ST_Buffer(ST_Buffer(way, 300, 'join=miter'), -300, 'join=miter')));

  1. Since I always had some “holes”, I filled them up with some PostGIS functions:

          (SELECT id, ST_Collect(ST_MakePolygon(way)) As way
           FROM (
               SELECT id, ST_ExteriorRing((ST_Dump(way)).geom) As way
               FROM city_boundaries
               WHERE (ST_Area(ST_Transform(way, 4326)::geography) / 1000000) > 3.5
               ) s
           GROUP BY id
          ) AS test

The result looks like very nice and I think it’s enough for my usage.

Thanks to all for your help!
Luca