This is the third post in a series about z-quads, a spatial coordinate system. The first one was about the basic construction, the second about how to determine which quads contain each other. This one is about conversions: how to convert positions to quads and back.

First, a quick reminder. A quad is a square(-ish) area on the surface of the earth. Most of the time we’ll use the unit square instead of the earth because it makes things simpler and you can go back and forth between lat/longs and unit coordinates just by scaling the value up and down. (There’s a section at the end about that.)

So, when I’m talking about encoding a point as a quad, what I mean more specifically is: given a point on the unit square (*x*, *y*) and a zoom level *z*, find the quad at zoom level *z* that contains the point. As it turns out, the zoom level doesn’t really matter so we’ll always find the most precise quad for a given point. If you need a different zoom level you can just have a step at the end that gets the precise quad’s ancestor at whatever level you want.

## Precision

So far we’ve focused on quads as an abstract mathematical concept and there’s been no reason to limit how far down you can zoom. In practice though you’ll want to store quads in integers of finite size, typically 32-bit or 64-bit. That limits the precision.

It turns out that you can fit 15 zoom levels into 32 bits and 31 into 64. At zoom level 15 the size of a quad is around 1.5 km^{2}. At zoom level 31 it’s around 3 cm^{2}. So yeah, you typically don’t need more than 31 levels. A third case is if you want to store your quad as a double – JavaScript, I’m looking at you here. A double can hold all integers up to 53 bits, after which some can’t be represented. That corresponds to 26 zoom levels which is also plenty precise for most uses.

## Encoding points

Okay, now we’re ready to look at encoding. What we’ll do is take a unit point (*x*, *y*) and find the level 31 quad that contains it. Since we’ve fixed the zoom level we’ll immediately switch to talking about the scalar rather than the quad. That makes everything simpler because within one level we’re dealing with a normal z-order curve and conversion to and from those are well understood.

The diagram on the right illustrates something we’ve seen before: that scalars are numbered such that the children at each level are contiguous. In this example, the first two green bits say which of the four 2×2 divisions we’re in. The next two say which 1×1 subdivision of that we’re in. It’s useful to think of this the same way as the order of significance within a normal integer: the further to the left a pair of bits is in a scalar the more significant it is in determining where the scalar is.

Now things will have to get really bit-fiddly for a moment. Let’s look at just a single 2×2 division. Each subdivision can be identified with an (*x*, *y*) coordinate where each ordinate is either 0 or 1, or alternatively with a scalar from 0 to 3. The way the scalars are numbered it happens that the least significant scalar bit gives you the *x*-coordinate and the most significant one gives you *y*:

(This is by design by the way, it’s why the z-shape for scalars is convenient). This means that given point *x* and *y*, each 0 or 1, you get the scalar they correspond to by doing

*x* + 2 *y*

What’s neat about that is that because of the way pairs of bits in the scalar increase in significance from right to left matches the significance of bits in an integer, you can do this for each bit of *x* and *y* at the same time, in parallel. To illustrate how this works I’ll take an example.

Let’s start with the point (2/5, 2/3) and find the level 5 scalar that contains it. First we’ll multiply both values by 3^{2} (=2^{5}) and round the result to an integer,

x = (2/5) × 32 ≈ 12 = 01100b y = (2/3) × 32 ≈ 21 = 10101b

Now we’ll perform *x* + 2 *y* for each bit in parallel like so,

Looks bizarre right, but it does work. We spread out the bits from each input coordinate so they can be matched up pairwise and then add them together. The result is the scalar we’re looking for because the order of significance of the bit-pairs in the scalar match the significance of the individual bits in each of the inputs.

Now, how can you tell that the result, scalar 626 on zoom level 5, is the right answer? This diagram suggests that it’s about right:

The bit-pairs in the scalar are [2, 1, 3, 0, 2] so we first divide and take subdivision 2, then divide that and take 1, then divide that, and so on. Ultimately we end up at the green dot. In this case the last step would be to add the appropriate bias, *b*_{5}, to get a proper quad and then we’re all done.

Because this conversion handles all the levels in parallel there’s little benefit to converting to a higher zoom level than the most precise, typically 31. If you want level 5 do the conversion to level 31 and then get the 26th ancestor of the result.

Most of the work in performing a conversion is spreading out the scaled 31-bit integer coordinates. One way to spread the bits is using parallel prefix which takes log_{2} *n* steps for an *n*-bit word,

00000000abcdefgh → 0000abcd0000efgh → 00ab00cd00ef00gh → 0a0b0c0d0e0f0g0h

That’s just one way to do it, there are others.

## Decoding quads

The operation above gets you into quad-land and obviously you’ll want to be able to get out again. This turns out to just be a matter of running the encoding algorithm backwards. The main difference is that here we’ll stay on whatever zoom level the quad’s on instead of always going to level 31.

Given a quad *q* this is how you’d decode it:

- Subtract the bias to get the quad’s scalar.
- Mask out all the even-index bits to get the spread
*x*value and the odd-index bits to get the spread*y*value. Shift the*y*value right by 1. - Pack the spread-
*x*and –*y*values back together to get a*z*-bit integer value for each. Packing can be similarly to spreading but in reverse._{q} - Floating-point divide the integer values by 2
. You now have your (^{zq}*x*,*y*) coordinates between 0 and 1.

A quad is an area but the result here is a point, the top left corner of the quad. This is often what you’re interested in, for instance if you’re calculating is the full area of the quad. There you can just add the width and height of the quad (which are both 2^{–zq}) to get the other corners. Alternatively you might be interested in the center of the quad instead. To get you can modify the last step of the algorithm to do,

Apply (2 *v* + 1) / 2^{zq+1} as a floating-point operation to each of the integer coordinates.

This is just a circuitous way to add 0.5 to *v* before dividing by 2* ^{zq}* but this way we can keep the value integer until the final division. If you want to minimize the loss of precision when round-tripping a point through a quad you should use the center since it’s the point with the smallest maximum distance to the points that belong to that quad.

So, now we can convert floating-point coordinates to quads and back again with ease. Well, relative ease. In the last part I’ll talk briefly about something I’ve been hand-waving until now: converting from lat/long to the unit square and back.

## Latitude and longitude

When you encounter a lat/long it is typically expressed using the WGS 84 standard. The latitude, a value between -90 and 90, gives you the north-south position. The longitude, a value between -180 and 180, gives you east-west.

A typical world map will look something like the one on the right here with America to the left and Asia to the right. The top left corner, northwest of North America, is lat 90 and long -180. The bottom right, south of New Zealand, has lat -90 and long 180. You could convert those ranges to an (*x*, *y*) on the unit square in any number of ways but the way I’ve been doing it follows the z-order curve the same way I’ve been using in this post: North America is in division 0, Asia in division 1, South America in division 2, and Australia in division 3. The conversion that gives you that is,

*x* = (180 + *long*) / 360

*y* = (90 – *lat*) / 180

This puts (0, 0) in the top left corner and (1, 1) in the bottom right and assigns the quads in that z-shaped pattern.

There’s a lot to be said about map projections at this point. I limit myself to saying that while it’s convenient that there’s a simple mapping between the space most commonly used for geographic coordinates and the space of quads, the Mercator projection really is a wasteful way to map a sphere onto a square. For instance, a full third of the space of quads is spent covering the Southern Ocean and Antarctica even though they take up much less than a third of the surface of the globe. I’d be interested in any alternative mappings that make more proportional use of the space of quads.

Okay, so now there’s a way to convert between WGS 84 and quads. I suspect there’s a lot more you can do with z-quads but this, together with the operations from the two previous posts, what I needed for my application so we’re at to the end of my series of posts about them. Well, except for one final thing: string representation. I’ve always been unhappy with how unmemorable geographic positions are so I made up a conversion from quads to sort-of memorable strings which I’ll describe in the last post.

In the meantime I made a little toy site that allows you to convert between positions and quads at a given zoom level. Here is an example,

This is quad 167159423 at zoom level 14 which covers, as always, Århus. The URL is

http://zquad.p5r.org?latlng=56.1676,10.2062&zoom=14

which kind of explains itself. You can also do

http://zquad.p5r.org?address=San+Francisco&zoom=14

Finally to see where a particular quad is try

http://zquad.p5r.org?quad=15386

If you want to explore the space of quads you can play around with the parameters yourself.

Pingback: Z-Quads: Containment | Hummus and Magnets

Pingback: Z-Quads: Strings | Hummus and Magnets

Alternatively, you can bitwise unpack floats to extract the mantissa & exponent directly then right shift for negative exponents if needed. Take the most-significant bits equivalent to your zoom and you get the bits (e.g., 01100 & 10101) directly. I think it’s more elegant than float*int including calculation of power of 2 with integer truncation (e.g., (2/5)*2^5=12). After all, this calculation you are doing *IS* exactly extracting the most-significant bits of the mantissa.

Right?

PS: I haven’t really compared this to what you do but I thought I’d throw it out there anyway.

That’s a good point, I hadn’t considered that. I would expect it to be pretty processor specific which way is fastest but certainly if you can bypass the float stack completely that would buy you a lot on some architectures.

One wrinkle that complicates this is that while shifting right by minus the exponent would just work up to 63 you have to be careful to saturate at 63, at least on intel, since shr.w just masks off the lowest 6 bits regardless of what the rest are. I’m sure there’s a nifty way to accomplish that, ideally without branching, but I haven’t looked into that.