cancel
Showing results for
Did you mean:

Team Exasol

# Background

The following is an implementation of the processes and formula described at Microsoft help page Bing Maps Tile System.

The functions can be used to

• Transform global GPS coordinates (WGS 84 / lat+long) into bing maps pixel coordinates
• Transform pixel coordinates into map tile coordinates
• Transform map tile coordinates into quadtree tile identifiers

We chose to implement all functions using builtin (free of extra charge) functionality:

• regular user-defined functions (aka pl/SQL) instead of Lua/Python/Java UDF.
• No GEOSPATIAL functions (ST_TRANSFORM etc.)

However, this forces us to either

• Use multiple functions when multiple values (X and Y coordinates) are required
• or use string-encoding where this happens.

We use the latter approach (formatting coordinates as "X/Y") for more compact code (less functions), and for easier migration towards geospatial data types.

# How to implement Bing maps and QuadTree calculation

## Step 1

### WGS-84 to Pixels

This is a straightforward implementation of the given Mercator-like projection onto a pixel grid of the given detail level (width/height = 256 * 2^detail).
Like in the original article, latitude, and longitude are clipped at ~85 (resp. 180) degrees in both directions.

```create function coords2pixel(latitude double, longitude double, detail int)
returns varchar(100)
is
sinLatitude double;
pixelX decimal(36,0);
pixelY decimal(36,0);
begin
latitude := greatest( -85.05112878, least( 85.05112878, latitude ) );
longitude := greatest(-180, least( 180, longitude ) );
sinLatitude := sin(latitude * pi() / 180);

pixelX := ((longitude + 180) / 360) * 256 * power(2,detail);
pixelY := (0.5 - log((1 + sinLatitude) / (1- sinLatitude)) / (4 * pi())) * 256 * power(2,detail);

return pixelX || '/' || pixelY;
end;
/
```

Please note that this will perform a virtual transposition of coordinate ordering: input is LAT/LONG (or vertical/horizontal when looking at a map) while the output is X/Y (horizontal/vertical).

## Step 2

### Pixels to Map Tile

As we know that each map tile has a size of 256x256 pixels, transforming pixel coordinates into tile coordinates is quite easy.
We use regular expressions with lookahead and look behind qualifiers for string parsing and divide each coordinate by 256:

```create function pixel2tile(pixel varchar(100))
returns varchar(100)
is
begin
return
floor( regexp_substr(pixel, '.*(?=/)') / 256 )
|| '/' ||
floor( regexp_substr(pixel, '(?<=/).*') / 256 );
end;
/
```

REGEXP_SUBSTR is used as a replacement for ST_X and ST_Y, which are part of the GEOSPATIAL license package.

## Step 3

### Map Tile to QuadTree Identifier

The overall process of interleaving coordinates using three different number systems sounds quite complicated, however, it boils down to:

• Interpret your input coordinates as binary data (examining single bits using BIT_CHECK)
• Putting these binary streams side-by-side, each cross-stream pair of bits creates one digit of your output coordinate/string.

This makes the conversion routine itself rather short:

```create or replace function tile2quad(tile varchar(100), detail int)
returns varchar(30)
is
xCoord int;
yCoord int;
lod int;
begin
xCoord := regexp_substr(tile, '.*(?=/)');
yCoord := regexp_substr(tile, '(?<=/).*');
for lod := 1 to detail do
end for;

end;
/
```

## Step 4

We now have a set of cascading functions that can transform GPS coordinates into Bing map tile identifiers. Let us check out Nuremberg (49.45N, 11.08E) at detail level 3:

```select coords2pixel(49.45, 11.08, 3);
--> '1087/699'

select pixel2tile('1087/699');
--> '4/2'

--> '120'
```

And indeed, according to the image on the link given above, Nuremberg should be inside mapped tile "120" (which is fully contained in the larger map tiles "12" and "1").

Short version for a smaller map tile with more details:

```select tile2quad(pixel2tile(coords2pixel(49.45, 11.08, 10)), 10);
--> '1202033313'
```