ST_SquareGrid for a Dynamic Street Map Index

OpenStreetMap + PostGIS + 0 Intermediate Layers

The Task

Recreate the Yorkville Streets Map.

At Kendall County, we offer GIS services to our local municipalities, and the City of Yorkville recently signed on with us. One of the tasks that came out of that was to create an updated version of their “streets map”, something the Police Department likes to have printed out.

The most recent copy was from 2022, and there have been some annexations and developments since then, so it was in need of refreshing.

Essential Components

• City Streets
• Municipal boundary
• An alphanumeric grid
• Street index based on grid
• Other nearby features for context

They sent along a copy of the most recent map, and it was a doozy. A 40MB PDF, which nearly all of the County’s parcels embedded as a vector layer, lots of other stuff.

But as long as the end result looks right, I can approach this my own way. So let’s look at this map and identify the essential things.

Being a street map, it obviously needs to include all the roads in the City. We want the city boundary shown here as well. The key feature of this map, though, is the alphanumeric grid. Like battleship! We want to look at the map, say “I’m looking for Spring Street”, check the “S” section, and get a list of grid cells we should look in.

Beyond that, we’ll want surrounding features for context, things like other cities, waterways, etc.

In video: show close-up pan of map.

Can I Do This with OpenStreetMap Data?

• Streets
• City boundary
• Index grid?

Whenever I start a new project, I always ask myself this question. Can I do this with OSM? If the features are missing, can we add them?

At the County, we use OSM for everything. We like to think we’ve got the place pretty well-mapped. Boundaries and streets are OSM’s bread and butter, and we actively work to keep those features current and high-quality.

But the grid… That’s not something that can be mapped in OSM, nor should it be.

Can I avoid creating any new layers?

• Map is being made in QGIS
• Already writing SQL query layers for the OSM data

There’s gotta be a function for this…

Another question I ask when I start a new project: can I get this done without creating new layers? I hate creating new layers, especially single-purpose ones.

This is our first time making this map for Yorkville, but odds are, the thing could get updated on an annual or as-needed basis. If the layers can all be based on OSM or some kind of spatial SQL, future maps would be that much simpler.

This was not something I’d done before, but I knew enough PostGIS to know that I could write an expression to just make ad-hoc geometry. Maybe I could use an envelope derived from the city shape and create some polygons?

The Grid

Step One: Research

Thanks, DuckDuckGo!

I did what any good GIS person does when faced with a novel problem: I did a web search!

I was expecting some kind of Stack Exchange answer or a Python tutorial, but I was utterly delighted to see the first search result was the PostGIS docs themselves!

Reading through the examples, it was exactly what I wanted. I could feed it the OSM Yorkville boundary, and it would spit out a grid, and in our local coordinate system!

Step Two: SQL

select (st_squaregrid(1000, ST_Transform(the_geom, 3435))).*
from kendall_areas
where tags ->> 'admin_level' = '8' and tags ->> 'name' = 'Yorkville'

Too many cells! We’ll have to tweak until we find the sweet spot of small enough grid cells to be useful, but big enough that the index doesn’t have too many values.

Side note: when I created this database, I called the shape column the_geom? I don’t know why, but it’s stuck since.

Step Three: Proper Grid Size

select (st_squaregrid(10000, ST_Transform(the_geom, 3435))).*
from kendall_areas
where tags ->> 'admin_level' = '8' and tags ->> 'name' = 'Yorkville'

Step Three: Proper Grid Size

select (st_squaregrid(5000, ST_Transform(the_geom, 3435))).*
from kendall_areas
where tags ->> 'admin_level' = '8' and tags ->> 'name' = 'Yorkville'

Step Three: Proper Grid Size

select (st_squaregrid(4800, ST_Transform(the_geom, 3435))).*
from kendall_areas
where tags ->> 'admin_level' = '8' and tags ->> 'name' = 'Yorkville'

Cell Labels

Good, but what are these i and j fields? I need letters and numbers! I’m sure there’s a perfectly good reason for them to start at 193 and 372.

Also, we want them sort of battleship-style, with A1 being at the top left.

We want 383 == 1, 382 == 2, etc., and 193 == A, 194 == B on the other columns.

BUT! I don’t want to hard-code anything. We could just stick a i - 193, but what if there’s an annexation, and the leftmost i value is 192? Or what if we change the grid size? We need everything to be dynamic.

Step Four: Proper Labels

with grid as (
  select (st_squaregrid(4800, ST_Transform(the_geom, 3435))).*
  from kendall_areas
  where tags ->> 'admin_level' = '8' and tags ->> 'name' = 'Yorkville'
),

indices as (
  select
  min(grid.i) as min_i,
  max(grid.j) as max_j
  from grid
)

select
  chr(i - indices.min_i + 65) as ltr,
  (j - indices.max_j - 1) * -1 as num,
  row_number() over() as id,
  geom
from grid
join indices on 1=1

To get that, we just select the min i and max j from our grid in a CTE. Then we can just use those values together with each cell’s i and j vales.

On the letter side, we use chr to convert the number to a character, and using the proper offset starts us at ‘A’. For letters, we just flop the axis and offset to get the number.

And I’m going to change their names to ltr and num.

And voila! We can take this same expression, and feed it a different municipality, or a future Yorkville that’s changed its shape, and it still holds.

The Streets

Centerlines

select
  osm_id,
  the_geom,
  tags ->> 'name' as name
from kendall_ways
where tags ->> 'highway' ~* '.*ary.*|motorway.*|trunk.*|residential.*|unclassified.*'

For our OSM database, we just dump everything, tags and all, into a few tables with jsonb columns. Getting the streets is easy enough with the ->> accessor.

We include a regex filter to limit the layer to main traffic roads and links.

Keep it City-Centric

select
  osm_id,
  the_geom,
  tags ->> 'name' as name
from kendall_ways
where tags ->> 'highway' ~* '.*ary.*|motorway.*|trunk.*|residential.*|unclassified.*'
  and tags ->> 'operator' ~* '.*Yorkville.*'

We don’t need all of these indexed. Normally, to ask for “roads in Yorkville”, we’d be looking at spatial overlay. But thanks to our former intern Adam, every road in the County has an operator tag! So we can filter solely by attribute. I’m pretty sure the tag only says “City of Yorkville”, but I am using regex just in case there’s a stray “Yorkville” or “Yorkville Highway Department”.

Keep it City-Centric

select
  osm_id,
  the_geom,
  tags ->> 'name' as name
from kendall_ways
where tags ->> 'operator' ~* '.*Yorkville.*|.*County.*|.*Illinois.*'
  and tags ->> 'highway' is not null and tags ->> 'name' is not null

But we want major state / county roads that happen to be in the City, too!

Oh, and as it happens, there are a couple of named service roads. So we’ll just filter for things that have any highway tag as well as a name.

Keep it City-Centric

select
  osm_id,
  the_geom,
  tags ->> 'name' as name
from kendall_ways
where tags ->> 'operator' ~* '.*Yorkville.*|.*County.*|.*Illinois.*'
  and tags ->> 'highway' is not null and tags ->> 'name' is not null or tags ->> 'route' = 'road'

And a state-operated road (Bridge Street) isn’t the same as the route (Illinois 47), but the person looking at the map might want either. We dump road routes into our database as single linestring features, but they’re not highway features, so we just amend the query to grab road routes, too.

And for roads that only belong to a route, sometimes the road doesn’t have is own name. So it turns out we need to filter those out, otherwise we get null entries in the index later.

Routes highlighted red for reference.

Which Cells are they In?

ST_Intersects

with idx as (
  -- our grid query from earlier
)

select 
  ST_Transform(the_geom, 3435) the_geom,
  tags ->> 'name' as name,
  ltr,
  num
from public.kendall_ways
join idx on ST_Intersects(ST_Transform(the_geom, 3435), geom)
where tags ->> 'operator' ~* '.*Yorkville.*|.*County.*|.*Illinois.*'
  and tags ->> 'highway' is not null and tags ->> 'name' is not null or tags ->> 'route' = 'road'

It’s PostGIS Day, so I feel like ST_Intersects has to show up.

So now we have our grid query from earlier. This time, we’ll bring it back as a CTE. Then we do a second CTE of the roads. To intersect properly, we need them in the same spatial reference.

Rookie Mistake

Rookie mistake, I’ve got a duplicate of, for instance, County Highway 13, for every cell of the grid I intersected with.

Aggregate!

with idx as (
  -- our grid query from earlier
),

roads as (
  select 
    ST_Union(ST_Transform(the_geom, 3435)) the_geom,
    tags ->> 'name' as name
  from public.kendall_ways
  where tags ->> 'operator' ~* '.*Yorkville.*|.*County.*|.*Illinois.*'
    and tags ->> 'highway' is not null and tags ->> 'name' is not null or tags ->> 'route' = 'road'
  group by tags ->> 'name'
)

select
  row_number() over() id,
  name,
  string_agg(ltr || '-' || num, ', ') grid_index
  st_union(the_geom) the_geom
from roads
join idx on st_intersects(the_geom, geom)
group by name

As it happens, I’ve got to aggregate twice, because OSM roads are split any time attributes change, so we end up having multiple Spring Streets, for instance, when the lane count changes.

Just to make things easier to keep clear, we’ll do the first step as its own CTE and aggregate our roads by name, so that every indexable road is intersecting the grid once.

Then we’ll do the second intersection separately. In that, we use string_agg to take the various index values and join them into a comma-delimited string.

Replacements

-- in the final select
  name label,
  case
    when name = 'West Street' then name
    when name like 'West %' then substring(name, 6) || ', West'
    when name like 'East %' then substring(name, 6) || ', East'
    when name like 'North %' then substring(name, 7) || ', North'
    when name like 'South %' then substring(name, 7) || ', South'
    else name
  end name,

The whole point of the index is to look up streets by name, not by prefix. But we do want them to be labeled with the prefix out front. A bit of find and replace into a new field for that.

There’s always an edge case, too. Have to hard-code “west street”.

Route Names

-- in the roads CTE
case
    when tags ->> 'name' is not null then tags ->> 'name'
    when tags ->> 'network' = 'US:US' then concat('US ', tags ->> 'ref')
    when tags ->> 'network' = 'US:IL' then concat('IL ', tags ->> 'ref')
  end as name

Back in the roads CTE, we need to do something about the names. For road routes, we don’t have names, so let’s use the network and ref tags to build a stand-in.

Other Stuff

That’s it for the big stuff, now we just need some visual context.

Municipalities

select
  osm_id,
  the_geom,
  tags ->> 'name' name
from kendall_areas
where tags ->> 'admin_level' = '8'

We want to highlight Yorkville, but also show the nearby municipalities. An easy OSM query, and some symbology.

We’ll duplicate the layer and make one of them an inverted polygon to sort of “fade out” things outside of the city.

Other Roads

A couple things here. One, we want to display all the roads, but restrict the index to Yorkville roads. Further, we want to label every city road, but only label major roads outside of the city.

To accomplish this, we do a couple things.

1. The grid-indexed roads are not sybmolized.
2. A non-index query layer of roads is added, and includes the operator and highway tags.
3. Symbology and labeling rules are made in QGIS using these fields.

The Map

Honestly, the query was the easy part. Getting the different pieces of the index lined up nicely was the most time-consuming part of all this. And then they asked for it to be portrait, not landscape.

The long and short of it is that I completely accomplished my task. I had a map of Yorkville’s roads, indexed the way they wanted, and mapped out accordingly.

But over and above that, I’d done it entirely with OpenStreetMap data and dynamic, query-driven layers.

Already since exporting this map, several annexations have come through, and new neighborhoods and roads are being built. When the request comes in to update this map, all I’ll need to do is open the project and re-export.

Well, and also move around all those index boxes to get them to fit nicely.

You’re welcome, future me.