Map of Birds Eye Coverage in EMEA & APAC (Part 3) [Archive]

Johannes Kebeck's Blog

09-15-2007, 03:25 PM

Step 2: Which tiles are covered by the shape at level 14? First we use SQL Server 2008 to determine the minimum bounding rectangle, for example:<pre>SELECT Shape.STEnvelope().STAsText() FROM FromBlom WHERE ID=24</pre><a href="http://11011.net/software/vspaste"></a>
In this example STEnvelope() calculated the minimum bounding rectangle and STAsText() provides the result as so called Well Known Text (WKT) which is one of the OGC-specifications we support. Besides that we support also Well Known Binary (WKB) and the Geographic Markup Language (GML).
The*WKT for such a polygon looks like this:
POLYGON ((52.150335311889648 0.017314910888671875, 52.309449434280396 0.017314910888671875, 52.309449434280396 0.25343132019042969, 52.150335311889648 0.25343132019042969, 52.150335311889648 0.017314910888671875))
<a href="http://blufiles.storage.msn.com/y1pLfDzppuxasgsapM5mZgx_bB9m-SjQycZnw3n6y_Z3Su6O7WsXf_ZCr5cN3XneEbi1k-8I4rgbOY"><img style="border-right:0px;border-top:0px;border-left:0px;border-bottom:0px" height=292 alt=image src="http://blufiles.storage.msn.com/y1pLfDzppuxasgtlpHZJtPyptI-tQPSiJ1niDB4pog_p9hGzDIyWPuO1_6YjNPQBu-wqQdG9EaIl-o" width=400 border=0></a>
Once we have the bounding rectangle we can start to determine the tiles in it:
<ol>
<li>We determine the locations of the nodes in the bounding rectangle<pre>SELECT Shape.STEnvelope().STPointN(1).STAsText() FROM FromBlom WHERE ID=24</pre>
STPointN(i) gives me a simple feature of type POINT like 'POINT (52.150335311889648 0.017314910888671875)' for a point of my choice. We will only need point 2 and 4 for our calculations. <a href="http://blufiles.storage.msn.com/y1pLfDzppuxasjYBUEqn0_Q7wI40AGMtUtBvDQb-LvYpGG-ZhOKOGJRo1NxGaIYvqJPfDmoszGdPNQ"><img style="border-right:0px;border-top:0px;border-left:0px;border-bottom:0px" height=412 alt=image src="http://blufiles.storage.msn.com/y1pLfDzppuxasikT1KfnrufajVykPbB5_zPoWvUF5b6dQKCZWT 8SghGusSxG-HijcneIS5Myqf2EaI" width=400 border=0></a> 
<li>Now we determine the latitude and longitude with the functions <pre>SELECT Shape.STEnvelope().STPointN(1).STX FROM FromBlom WHERE ID=24
SELECT Shape.STEnvelope().STPointN(1).STY FROM FromBlom WHERE ID=24 </pre><pre>The result is of data type float e.g. 52.1503353118896 and 0.0173149108886719</pre>
<li>
Well and then we use some mathematical functions which are well described in Joe Schwarz's <a href="http://msdn2.microsoft.com/en-us/library/bb259689.aspx">article on MSDN</a> to determine the quadkey of these tiles. I implemented these functions in the following T-SQL script:<pre>DECLARE @ID INT;
DECLARE @i INT;
DECLARE @Lat FLOAT;
DECLARE @Lon FLOAT;
DECLARE @lvl INT;
DECLARE @pixelX FLOAT;
DECLARE @pixelY FLOAT;
DECLARE @tileXdecMin INT;
DECLARE @tileYdecMin INT;
DECLARE @tileXdecMax INT;
DECLARE @tileYdecMax INT;
DECLARE @tileXdec INT;
DECLARE @tileYdec INT;
DECLARE @tileXdecDummy INT;
DECLARE @tileYdecDummy INT;
DECLARE @tileXbin VARCHAR(38);
DECLARE @tileYbin VARCHAR(38);
DECLARE @j INT;
DECLARE @quadKey VARCHAR(19);

SET @lvl = 14;

DECLARE IdCursor CURSOR FOR SELECT ID FROM FromBlom WHERE Country='IE' ORDER BY ID;
OPEN IdCursor;
FETCH NEXT FROM IdCursor INTO @ID;
WHILE @@FETCH_STATUS = 0
BEGIN
 SET @i = 2;
 WHILE @i <= 4 BEGIN
 -- Get latitude and longitude of the bounding reactangle
 SET @Lat = (SELECT Shape.STEnvelope().STPointN(@i).STX FROM FromBlom WHERE ID=@ID);
 SET @Lon = (SELECT Shape.STEnvelope().STPointN(@i).STY FROM FromBlom WHERE <a href="mailto:ID=@ID);">ID=@ID);
</a> -- Calculate pixel X and Y of the CP
 SET @pixelX = (@Lon + 180) / 360 * 256 * POWER(2, @lvl);
 SET @pixelY = (0.5 - (LOG((1 + SIN(@Lat * PI() / 180)) / (1 - SIN(@Lat * PI() / 180))) / 4 / PI())) * 256 * POWER(2, @lvl);
 -- Calculate decimal tile X and Y
 IF (@i=2) BEGIN
 SET @tileXdecMin = FLOOR(@pixelX / 256);
 SET @tileYdecMin = FLOOR(@pixelY / 256);
 END;
 ELSE BEGIN
 SET @tileXdecMax = FLOOR(@pixelX / 256);
 SET @tileYdecMax = FLOOR(@pixelY / 256);
 END;
 SET @i = @i+2;
 END;

-- Now we loop through the tiles in this rectangle
 SET @tileYdec = @tileYdecMin;
 WHILE @tileYdec <= @tileYdecMax BEGIN
 SET @tileXdec = @tileXdecMin;
 WHILE @tileXdec <= @tileXdecMax BEGIN
 SET @tileYdecDummy = @tileYdec;
 SET @tileXdecDummy = @tileXdec;
 -- Calculate binary tile X and Y
 SET @tileXbin = '';
 SET @tileYbin = '';
 SET @j = @lvl-1;
 WHILE @j >= 0 BEGIN
 SET @tileXdecDummy = @tileXdecDummy - POWER(2, @j);
 IF (@tileXdecDummy < 0) BEGIN
 SET @tileXbin = @tileXbin + '0'
 SET @tileXdecDummy = @tileXdecDummy + POWER(2, @j);
 END
 ELSE BEGIN
 SET @tileXbin = @tileXbin + '1'
 END

 SET @tileYdecDummy = @tileYdecDummy - POWER(2, @j);
 IF (@tileYdecDummy < 0) BEGIN
 SET @tileYbin = @tileYbin + '0'
 SET @tileYdecDummy = @tileYdecDummy + POWER(2, @j);
 END
 ELSE BEGIN
 SET @tileYbin = @tileYbin + '1'
 END
 SET @j = @j-1;
 END;
 -- Calculate Quadkey
 SET @j = 1;
 SET @quadKey = '';
 WHILE @j <= @lvl BEGIN
 SET @quadKey = @quadKey + CAST(2 * CAST(SUBSTRING(@tileYbin, @j, 1) AS INT) + 1 * CAST(SUBSTRING(@tileXbin, @j, 1) AS INT) AS VARCHAR(1));
 SET @j = @j+1;
 END;
 INSERT INTO [VE-Tiles] (ID, LVL, Quadkey) VALUES (@ID, @lvl, @quadKey)
 SET @tileXdec = @tileXdec+1;
 END;
 SET @tileYdec = @tileYdec+1;
 END;
 PRINT '';
FETCH NEXT FROM IdCursor INTO @ID;
 END;
CLOSE IdCursor;
DEALLOCATE IdCursor;
<a href="http://blufiles.storage.msn.com/y1pLfDzppuxasjrQHjIC0LWIrc49iiCXKAC4EEgTT4T514o7U7 Uo2WPNr13A0mBVUZO2xbiODXkReU"><img style="border-right:0px;border-top:0px;border-left:0px;border-bottom:0px" height=340 alt=image src="http://blufiles.storage.msn.com/y1pLfDzppuxasiIgOg4RjvUKH0f0ewUMlmKEM36X4AxANGkV2K C-cjoW8djwmOx7FmhJzd1KhVHkrM" width=400 border=0></a> </pre>
If you had a closer look at my script you noticed that I inserted the quadkeys in a new table and in the next step I will create a shape which represents the bounding box of each tile:<pre>DECLARE @myQuadkey VARCHAR(19);
DECLARE @myDummy VARCHAR(38);
DECLARE @myLevel INT;
DECLARE @myQuadKeyBin VARCHAR(38);
DECLARE @i INT;
DECLARE @TileXBin VARCHAR(19);
DECLARE @TileYBin VARCHAR(19);
DECLARE @TileX INT;
DECLARE @TileY INT;
DECLARE @PixelXMin INT;
DECLARE @PixelXMax INT;
DECLARE @PixelYMin INT;
DECLARE @PixelYMax INT;
DECLARE @LatMin FLOAT;
DECLARE @LatMax FLOAT;
DECLARE @LongMin FLOAT;
DECLARE @LongMax FLOAT;
DECLARE Tile_Cursor CURSOR FOR SELECT Quadkey FROM [VE-Tiles] WITH (NOLOCK) WHERE BBox IS NULL;
OPEN Tile_Cursor;
FETCH NEXT FROM Tile_Cursor INTO @myQuadkey;
WHILE @@FETCH_STATUS = 0
BEGIN
-------------------------------------------------------------------
-- Convert Quadkey in binary Quadkey
-------------------------------------------------------------------
 SET @myLevel = LEN(@myQuadkey);
 SET @myQuadKeyBin = '';
 SET @i = 0;
 WHILE @i < @myLevel BEGIN
 SET @myDummy = SUBSTRING(@myQuadkey, @myLevel-@i, 1);
 IF @myDummy = 3 SET @myQuadKeyBin = '11' + @myQuadKeyBin;
 IF @myDummy = 2 SET @myQuadKeyBin = '10' + @myQuadKeyBin;
 IF @myDummy = 1 SET @myQuadKeyBin = '01' + @myQuadKeyBin;
 IF @myDummy = 0 SET @myQuadKeyBin = '00' + @myQuadKeyBin;
 SET @i = @i+1
END
-------------------------------------------------------------------
-- Break binary Quadtree in binary tileX and tileY
-------------------------------------------------------------------
 SET @myDummy = @myQuadKeyBin;
 SET @TileXBin = '';
 SET @TileYBin = '';
 SET @i = 0;
 WHILE @i < LEN(@myQuadKeyBin) BEGIN
 SET @myDummy = SUBSTRING(@myQuadKeyBin, LEN(@myQuadKeyBin)-@i, 1)
 IF (CAST(@i AS FLOAT)/2 - FLOOR(@i/2) = 0) BEGIN
 SET @TileXBin = @myDummy + @TileXBin;
 END
 ELSE BEGIN
 SET @TileYBin = @myDummy + @TileYBin;
 END
 SET @i = @i+1
END
-------------------------------------------------------------------
-- Convert binary tileX and tileY to decimal tileX and tileY
-------------------------------------------------------------------
 SET @TileX = '';
 SET @TileY = '';
 SET @i = 0;
 WHILE @i < @myLevel BEGIN
 SET @myDummy = SUBSTRING(@TileXBin, LEN(@TileXBin)-@i, 1)
 SET @TileX = @TileX + @myDummy * POWER(2, @i)
 SET @i = @i+1
END
SET @i = 0;
 WHILE @i < @myLevel BEGIN
 SET @myDummy = SUBSTRING(@TileYBin, LEN(@TileYBin)-@i, 1)
 SET @TileY = @TileY + @myDummy * POWER(2, @i)
 SET @i = @i+1
END
-------------------------------------------------------------------
-- Determine PixelX and PixelY-Ranges
-------------------------------------------------------------------
 SET @PixelXMin = @TileX * 256
SET @PixelXMax = ((@TileX + 1) * 256) - 1
SET @PixelYMin = @TileY * 256
SET @PixelYMax = ((@TileY + 1) * 256) - 1
-------------------------------------------------------------------
-- Calculate Bounding-Box
-------------------------------------------------------------------
 SET @LongMin = (CAST(@PixelXMin AS FLOAT) / 256 / POWER(2, CAST(@myLevel AS FLOAT)) * 360) - 180
SET @LongMax = (CAST(@PixelXMax AS FLOAT) / 256 / POWER(2, CAST(@myLevel AS FLOAT)) * 360) - 180
SET @LatMin = ASIN((EXP((0.5 - CAST(@PixelYMin AS FLOAT) / 256 / POWER(2, CAST(@myLevel AS FLOAT))) * 4 * PI()) - 1) / (EXP((0.5 - CAST(@PixelYMin AS FLOAT) / 256 / POWER(2, CAST(@myLevel AS FLOAT))) * 4 * PI()) + 1)) * 180 / PI()
 SET @LatMax = ASIN((EXP((0.5 - CAST(@PixelYMax AS FLOAT) / 256 / POWER(2, CAST(@myLevel AS FLOAT))) * 4 * PI()) - 1) / (EXP((0.5 - CAST(@PixelYMax AS FLOAT) / 256 / POWER(2, CAST(@myLevel AS FLOAT))) * 4 * PI()) + 1)) * 180 / PI()
-------------------------------------------------------------------
-- Update Geometry
-------------------------------------------------------------------
 EXEC('UPDATE [VE-Tiles] SET BBox=geometry::STGeomFromText(' + '''' + 'POLYGON((' + @LatMin + ' ' + @LongMin + ', ' + @LatMin + ' ' + @LongMax + ', ' + @LatMax + ' ' + @LongMax + ', ' + @LatMax + ' ' + @LongMin + ', ' + @LatMin + ' ' + @LongMin + '))' + '''' + ', 4326) WHERE QuadKey=' + '''' + @myQuadkey + '''');
 FETCH NEXT FROM Tile_Cursor INTO @myQuadkey;
 END;
CLOSE Tile_Cursor;
DEALLOCATE Tile_Cursor;</pre><a href="http://11011.net/software/vspaste"></a><pre>With a call like the one shown below we can retrieve the result for one shape</pre><pre>SELECT BBox.STAsText() FROM [VE-Tiles] WHERE ID=24</pre><pre>and here is how it looks like in Virtual Earth:</pre><pre><a href="http://blufiles.storage.msn.com/y1pLfDzppuxasg2J8I4a7W9Lb5nZCdYNoH3uH2arMiWr7ZgQD-kw1OMTYFQIlQ6Cy-AEnYpDpTmg8I"><img style="border-right:0px;border-top:0px;border-left:0px;border-bottom:0px" height=291 alt=image src="http://blufiles.storage.msn.com/y1pLfDzppuxasiOHgKJg79GJ-ErGxoo8X7Cz6n7JLA82jGOVkr68UDOv5B9de2dnNiKmGulChro hmI" width=400 border=0></a> </pre></ol><img src="http://c.services.spaces.live.com/CollectionWebService/c.gif?cid=4819404664324524694&page=RSS%3a+Map+of+Birds+Eye+Coverage+in+EMEA+%26+ APAC+(Part+3)&referrer=" width="1px" height="1px" border="0" alt=""><img style="position:absolute" alt="" width="0px" height="0px" src="http://c.live.com/c.gif?NC=31263&NA=1149&PI=88469&RF=&DI=3919&PS=85545&TP=johanneskebeck.spaces.live.com&GT1=JohannesKebeck">

Click here to view the article. (http://JohannesKebeck.spaces.live.com/Blog/cns!42E1F70205EC8A96!1736.entry)