The SqlGeometry with Microsoft SQL Server

I came across a curious error earlier today when attempting to use a SqlDataReader to read a column of geometry data from a SQL Server table:

System.InvalidCastException: Unable to cast object of type ‘Microsoft.SqlServer.Types.SqlGeometry’ to type ‘Microsoft.SqlServer.Types.SqlGeometry’

SqlGeometry to SqlGeometry… you’d think that would be a pretty easy cast, wouldn’t you? It turns out that this is a problem caused by a conflict between the spatial libraries used in SQL Server Denali compared to that in 2008/R2, and you’ll get this error depending on which version of Microsoft.SqlServer.Types.dll you use, and how you try to access geometry or geography columns from a datareader:

[php]
while (dataReader.Read())
{
// This works in SQL Server 2008/R2, but errors with Denali
SqlGeometry g = (SqlGeometry)dataReader.GetValue(0);

// This works in SQL Server 2008/R2, but errors with Denali
SqlGeometry g = (SqlGeometry)dataReader["GeomCol"];

// This works in Denali, but not in SQL Server 2008/R2
SqlGeometry g = SqlGeometry.Deserialize(reader.GetSqlBytes(0));

// This works in Sql Server 2008/R2/Denali
SqlGeometry g = new SqlGeometry();
g.Read(new BinaryReader(reader.GetSqlBytes(0).Stream));
}
[/php]

After a bit of digging around, it appears that using GetValue or square brackets notation [] to access a geometry/geography field in a SqlDataReader is hard-coded to load the 10.0 (SQL Server 2008) version of the Microsoft.SqlServer.Types library.

If you’ve got side-by-side installations of both SQL Server 2008/R2 and Denali (as I have), and try to reference the 11.0 (Denali) version of Microsoft.SqlServer.Types, you’ll therefore get an assembly mismatch when both versions of the library are loaded, which causes the slightly unhelpful error listed at the top of this post. Even if you’ve only got Denali installed, your code may still try to reference a (non-existent) 2008/R2 version of the Microsoft.SqlServer.Types.dll library, so you’ll get a different error instead:

Could not load file or assembly ‘Microsoft.SqlServer.Types, Version=10.0.0.0, Culture=neutral, PublicKeyToken=89845dcd8080cc91′ or one of its dependencies. The system cannot find the file specified.

The simplest way to resolve these errors is by changing the way you reference any geography/geometry (and, I imagine, hierarchyid) columns from your DataReader, as in the code example above. Alternatively, you can set up an assembly redirection in the application configuration file as explained here (about halfway down), which will allow you to correctly target the Denali version.

As per the What’s new in SQL Server Denali whitepaper, “… side-by-side installations of SQL Server Code-Named “Denali” CTP1 and CTP3 are not supported with existing SQL Server 2008 installations …”, so perhaps I only have myself to blame for this. Interestingly though, the person who raised this MS Connect issue, says that they have experienced exactly the same problem on a clean install of Denali CTP3. The response from Microsoft suggests that this may be due to older versions of the library being packaged with Visual Studio 2010, and also confirms that the problem will not be resolved prior to RTM of SQL Server Denali.

Strangely, I encountered another curious error a few months ago concerning version conflicts of Microsoft.SqlServer.Types. My CTP3 Management Studio Spatial Results tab does not plot curved geometries (selecting a CircularString or the result of BufferWithCurves etc. just produces a blank pane). I had originally assumed that, since this was only a CTP release, this feature had simply not been added yet. It turns out that curved geometries are supported in SSMS CTP3 Spatial Results tab but, if you have side-by-side SQL Server 2008 and Denali, this can corrupt this feature. I guess the reason is similar – that SSMS is somehow attempting to load the SQL Server 2008/R2 version of Microsoft.SqlServer.Types, which, of course, doesn’t support curved geometries.

Deleting Blank Tiles

 

Creating raster tilesets almost invariably leads to the creation of some blank tiles – covering those areas of space where no features were present in the underlying dataset. Depending on the image format you use for your tiles, and the method you used to create them, those “blank” tiles may be pure white, or some other solid colour, or they may have an alpha channel set to be fully transparent.

Here’s an example of a directory of tiles I just created. In this particular z/x directory, more than half the tiles are blank. Windows explorer shows them as black but that’s because it doesn’t process the alpha channel correctly. They are actually all 256px x 256px PNG images, filled with ARGB (0, 0, 0, 0):

image

What to do with these tiles? Well, there’s two schools of thought:

  • The first is that they should be retained. They are, after all, valid image files that can be retrieved and overlaid on the map. Although they aren’t visually perceptible, the very presence of the file demonstrates that the dataset was tested at this location, and confirms that no features exist there. This provides important metadata about the dataset in itself, and confirms the tile generation process was complete. The blank images themselves are generally small, and so storage is not generally an issue.
  • The opposing school of thought is that they should be deleted. It makes no sense to keep multiple copies of exactly the same, blank tile. If a request is received for a tile that is not present in the dataset, the assumption can be made that it contains no data, and a single, generic blank tile can be returned in all such instances – there is no benefit of returning the specific blank tile associated with that tile request. This not only reduces disk space on the tile server itself, but the client needs only cache a single blank tile that can be re-used in all cases where no data is present.

I can see arguments in favour of both sides. But, for my current project, disk and cache space is at a premium, so I decided I wanted to delete any blank tiles from my dataset. To determine which files were blank, I initially thought of testing the filesize of the image. However, even though I knew that every tile was of a fixed dimension (256px x 256px), an empty tile can still vary in filesize according to the compression algorithm used. Then I thought I could loop through each pixel in the image and use GetPixel() to retrieve the data to see whether the entire image was the same colour, but it turns out that GetPixel() is slooooowwwww….

The best solution I’ve found is to use an unsafe method, BitMap.LockBits to provide direct access to the pixel byte data of the image, and then read and compare the byte values directly. In my case, my image tiles are 32bit PNG files, which use 4 bytes per pixel (BGRA), and my “blank” tiles are completely transparent (Alpha channel = 0). Therefore, in my case I used the following function, which returns true if all the pixels in the image are completely transparent, or false otherwise:

public static Boolean IsEmpty(string imageFileName)
{
  using (Bitmap b = ReadFileAsImage(imageFileName))
  {
    System.Drawing.Imaging.BitmapData bmData = b.LockBits(new Rectangle(0, 0, b.Width, b.Height), System.Drawing.Imaging.ImageLockMode.ReadOnly, b.PixelFormat);
    unsafe
    {
      int PixelSize = 4; // Assume 4Bytes per pixel ARGB
      for (int y = 0; y < b.Height; y++)
      {
        byte* p = (byte*)bmData.Scan0 + (y * bmData.Stride);
        for (int x = 0; x < b.Width; x++)         {           byte blue = p[x * PixelSize]; // Blue value. Just in case needed later           byte green = p[x * PixelSize + 1]; // Green. Ditto           byte red = p[x * PixelSize + 2]; // Red. Ditto           byte alpha = p[x * PixelSize + 3];           if (alpha > 0) return false;
        }
      }
    }
    b.UnlockBits(bmData);
  }

  return true;
}

 

It needs to be compiled with the /unsafe option (well, it did say in the title that this post was dangerous!). Then, I just walked through the directory structure of my tile images, passing each file into this function and deleting those where IsEmpty() returned true.

 

The DEM Shaded Overlays

 

The default Bing Maps road style uses a “hillshade” effect to give an impression of underlying terrain. It’s a relatively subtle, but surprisingly powerful technique to enhance the appearance of map layers, as demonstrated by comparing the following two tiles:

image
Without hillshading		 image
With hillshading

In this post, I’ll describe how to create your own hillshade overlay from digital elevation model (DEM) data, using the GDAL toolset.

By creating the overlay as a set of semi-transparent tiles, rather than pre-rendered into the tiles as shown above, you can place them on top of any Bing Maps/Google Maps et al. tilelayer to represent the underlying terrain.

The process I’ve followed is based on the work of others, most notably PerryGeo, and you can find some other guides on the internet to achieve this same effect. However, I found some of the existing guides on the subject to be either out-of-date or require knowledge of Linux BASH commands etc., so I hope that some of you will find this new step-by-step guide helpful.

1.) Acquire a DEM terrain model

To start with, you’re going to need some source data about the underlying terrain of the earth from which to calculate your hillshade. There’s lots of places to acquire this data from; Perhaps the easiest to use (assuming you’ve got Google Earth installed) is to open the kmz file available from http://www.ambiotek.com/topoview. This uses Google Earth as a graphical interface for v4.1 of the  elevation dataset gathered by the Shuttle Radar Topography Mission (SRTM), from which you can click to download individual DEM tiles covering 5°x 5°, as shown below:

Alternatively, you can access these files directly from the KCL server (my former university, incidentally) at http://srtm.geog.kcl.ac.uk/portal/srtm41/srtm_data_geotiff/

The data is provided in GeoTIFF format. You can load one of these tiles up in any graphics program that can load TIFF files, but it won’t look very interesting yet. The height information is encoded in additional metadata that will be ignored by normal graphics programs, so you’ll probably just get an image like this (this is srtm_36_02.tif):

Black parts show the presence of data in the underlying file, which we’ll subsequently process using GDAL tools to create shaded images.

2.) Reproject to Spherical Mercator

Most DEM data sources, including the SRTM data I linked to above, are provided in Plate Carree projection – i.e. WGS84 coordinates of longitude are mapped directly to the x axis of the image, while latitude is mapped directly to the y axis. Before we create tiles from this data suitable for overlay on Bing Maps, Google Maps, et al. we therefore need to transform it into the Spherical Mercator projection. You can do this using gdalwarp, as follows:











gdalwarp -dstnodata 0 -tr 305.7481 305.7481 -multi -co "TILED=YES" -t_srs EPSG:3857 srtm_36_02.tif srtm_36_02_warped.tif

The full list of parameters accepted by gdalwarp are listed here,  but the options I set are as follows:

  • dstnodata states what value to use to represent nodata values (the equivalent of null in a SQL database, for example). I’ve set a value of 0 (i.e. black).
  • tr gives the target resolution in x and y dimensions. The SRTM data I’m using was recorded at 90m resolution, so you might think that this should be set to 90 90. However, I’m going to be using this data for display on Bing Maps at different zoom levels, which will necessarily involve resampling the image.  Therefore, you should set this value to the resolution (in metres/pixel) of the maximum zoom level on which you plan to overlay your data. (Remember that maximum zoom level will have the smallest resolution). You can obtain this value from my Bing Maps Ready Reckoner. In the case above, I’m planning overlaying my data on Zoom Level 9 and above, so I set a value of 305.7481 (in both dimensions). If I’d wanted to go to Zoom Level 10, I would have decreased this to 152.87 instead.
  • multi allows parallel processing
  • co “TILED=YES” is a format-specific option that states that the output TIFF file should be tiled rather than stripped (see http://www.fileformat.info/format/tiff/egff.htm for an explanation of the difference)
  • t_srs gives the destination spatial reference system into which the image should be reprojected. In this case, EPSG:3857, as used by Bing Maps, Google Maps etc.

The resulting image, srtm_36_02_warped.tif, will still be a GeoTIFF file, but will now be projected as follows. The height and width of the output image will depend on the target resolution you specified in the tr parameter:

 

3.) Convert from DEM to Hillshade

The warped GeoTIFF file has height data encoded in it, but we want to translate that information into a visible shaded effect, and for this we can use gdaldem.

gdaldem actually provides several interesting functions related to working with DEM data, including the ability to derive contour lines, and create shaded relief maps. Maybe I’ll write about these another time, but for this example we want to use the hillshade mode. You can shade the warped image created in the previous step as follows:











gdaldem hillshade srtm_36_02_warped.tif srtm_36_02_warped_hillshade.tif -z 2 -co "TFW=YES"

This time, I’m only supplying two additional parameters:

  • z is a scaling factor applied to the generated hillside image that accentuates the hills, increasing the contrast of the image. I provided a value of 2 just to enhance the effect a bit, but you might decide you don’t need this.
  • co “TFW=YES” specifies that the output image should be created with an accompanying “world file”. This is a simple ASCII text file that provides additional information about the geographic extents of the created image, which we’ll need to use in a later step to line the hillshade image up with the Bing Maps tiling system. You can look up more information about world files on wikipedia.

There are additional parameters that allow you to specify the direction and the angle of the light source from which the simulated shadows will be created.

The result of executing the above code will be another TIFF file, in which the background is black, and the elevation data from the DEM has been converted into shades of grey, as follows:

 

At this stage, you could stop if you wanted to, and simply create a tile layer from the hillshaded image, which would look a bit like this:

 

Which makes the landscape of North Wales look a bit like the Moon, I think…

To make the data slightly more usable, we need to carry on with a few more tweaks.

4.) Making a Semi-Transparent Overlay

Currently, our hillshade image is opaque, with the shadows cast by terrain represented by variations in brightness of the colour used. To make this into an re-usable overlay that can be used on top of other layers, we need to make the image semitransparent, with shadows cast by terrain being represented by variations in opacity instead.

There are several ways of modifying the image data to achieve this effect. You could do it in Photoshop or another graphics program, for example, or using the graphics libraries in C# or PHP. Since I’m currently trying to learn Python, and GDAL is quite closely linked with Python, I’ll try to do it using the Python Imaging Library instead.

The following Python script makes a number of tweaks to the image above. Firstly, it converts it to a pure greyscale image (while the image above looks greyscale, it’s actually using a colour palette). It then inverts the image, turning it into a negative image. The reason for the inversion is that we then copy the (single) channel of the greyscale image into the opacity channel of a new RGBA image – areas that were very light in the source want to have very low opacity in the transparent image, and vice-versa, so the channel needs to be inverted.

Finally, we scan through the data to find instances of pixels that are pure black (RGBA value 0, 0, 0, 255) –this was the nodata value we set in step one – and replace them with pure transparent pixels (0, 0, 0, 0). The alpha channel in the tuples of any other pixels is also lightened slightly – I chose a value of 74 somewhat arbitrarily because I thought the resulting image looked good – you can choose whatever value you want, or none at all.
[php]
from PIL import Image as PImage
from PIL import ImageOps

# Load the source file
src = PImage.open("srtm_36_02_warped_hillshade.tif")

# Convert to single channel
grey = ImageOps.grayscale(src)

# Make negative image
neg = ImageOps.invert(grey)

# Split channels
bands = neg.split()

# Create a new (black) image
black = PImage.new(‘RGBA’, src.size)

# Copy inverted source into alpha channel of black image
black.putalpha(bands[0])

# Return a pixel access object that can be used to read and modify pixels
pixdata = black.load()

# Loop through image data
for y in xrange(black.size[1]):
for x in xrange(black.size[0]):
# Replace black pixels with pure transparent
if pixdata[x, y] == (0, 0, 0, 255):
pixdata[x, y] = (0, 0, 0, 0)
# Lighten pixels slightly
else:
a = pixdata[x, y]
pixdata[x, y] = a[:-1] + (a[-1]-74,)

# Save as PNG
black.save("srtm_36_02_warped_hillshade_alpha.png", "png")
[/php]
(Much of the logic in this script came from here). The resulting image will be a PNG file, in which darker shadows are represented by increasingly opaque black parts, while lighter shadows are more transparent: