SharpGIS

#GIS from a .NET developer's perspective

SharpMap v0.9 RC1 online !

SharpMap v0.9RC1 is now available for download from the SharpMap website.

Thanks to all the people who kept sending feedback and bug reports!

I've resolved a lot of bugs since beta 1, as well as added some obviously missing properties and methods in the API. Especially the Shapefile reader has been enhanced to handle a lot of special cases (in one case it was even able to render a file that ArcGIS claimed to be corrupt). On-the-fly transformation has also been completed and now supports transformation between geographic and projected coordinate systems.

Thanks to Humberto Ferreira for providing us with an Oracle Spatial provider as well. SharpMap now supports data from ESRI Shapefiles, PostGIS, Oracle, Microsoft SQL Server, Access and ECW raster.

Last but not least, I've spent a great deal of time enhancing the documentation, so go take a look at that as well. There are lots of code-samples hidden in there as well.

Make Microsoft SQL Server geospatial

I’ve always thought that on the spatial support, MSSQL was way behind many of the other database servers, in its lack of supporting storage of geometry. With the new .NET CLR you can actually add your own .NET-based object types and I’ve also tried implementing the Simple Features Specification in SQL Server. There are some limitations that made me give this up though. First of all, a user data type cannot be more than 8000 bytes. That is at most no more than 500 vertices in a geometry object, which is far too little for an accurate coastline for instance. Another problem is that SQL Server doesn’t support inheritance chains, so you can’t make a good object-oriented implementation of your datatype either.

…so yesterday I went for a completely different and much simpler approach. I decided to just store the geometry as Well-Known Binary in an image column. The reason for using an image column is that it can hold up to 2Gb of data, which should be sufficient for most geometry objects ;-). A binary field has the same 8000 byte limitation as UDT so this is no good. In addition to the geometry field, I create four real-type fields, holding the min/max values of the geometry envelope. This makes it efficient to do boundingbox based queries on the data. Any additional fields would be attributes for the geometry object.

I implemented the whole thing using SharpMap. First I created a small upload application that takes a shapefile, creates a table in the database and uploads the geometry and attributes to it. SharpMap has the necessary data readers and WKB formatters for this. The second part was to create a data provider from which SharpMap could draw from. I more or less based this on the PostGreSQL/PostGIS data provider for SharpMap, by changing the boundingbox query to use the four envelope fields. All this was not much more than an hour’s work, so it is very simple to accomplish.

I must say I was very surprised by the performance of the approach. It is just slightly faster than the shapefile data provider, which used to be the fastest data provider for SharpMap. In comparison the PostGreSQL/PostGIS is generally 4-6 times slower.

I have created a small demo web-application you can download from here. It contains two pages: one for uploading to the database, and one for rendering a layer from the database. All you need to do is to add an empty SQL Server 2005 Express database to the \App_Data\ folder and name id "GeoDatabase.mdf".

Download SharpMapSqlServer.zip (181,74 KB) (updated May 20, 2006)

Update: The MsSqlProvider is now also included in v0.9RC1, including a method for uploading from any of the SharpMap datasources to MS SQL.

Transformation problems with NetTopologySuite and GeoTools.NET

The last couple of days I've been working on doing on-the-fly transformations in SharpMap. I started out by implementing all the interfaces in the OGC "Coordinate Transformation Services" specification. These specifications are just great and makes my work of designing APIs so much eaier. Anyway when I started at the actual implementation I also looked at NetTopologySuite (which is a .NET port of Java Topology Suite and GeoTools.NET for a little help on the more complex transformation formulas.

It didn't really work out. All my results was either wrong, resulted in Overflow and DivisionByZero exceptions and NaN values. All I could then do was working through the EPSG transformation formulas, which was I was trying to avoid in the first place (oh by the way, EPSG recently released v7.10 of their coordinate system database).

Well I finally hacked it, and my conclusion is that there are several problems in the NTS code that causes very inaccurate (if not completely wrong) transformations. Problems like an incorrect constant for converting between radians and degrees, transformation formulas that are wrong (at least the Mercator projection - haven't checked the others yet), transformation-factories that doesn't respect the order of source and target coordinate systems (it will always transform from geographic to projected coordsys), and using the TransformList(array-of-points) method in many cases gives a different result than Transform(point).

I emailed my findings to Diego and Andrew who's doing a great job on NTS, so that they can benefit from this work well, but this is just to warn you people out there that NTS is probably giving your incorrect or inaccurate transformation results, until the bug is corrected. This is by no means Diego's fault; I guess he just copied it from GeoTools.NET (since the code is more or less the same), but this reminds me why we always use the national cadastras transformation libraries for transforming data...

Oh well as you can probably understand from this post, SharpMap now supports on-the-fly transforrmation! How cool is that? :-) If you can't wait for the next release, grap the source from the GotDotNet workspace. Warning though: I haven't checked the rest of the formulas thoroughly yet, but Mercator projection is checked and corrected and UTM seems to be working...

Link between NetTopologySuite and SharpMap

One of the things that SharpMap at the current stage lacks most is probably a relation model, making for instance intersectiontests possible (unless you use the PostGIS provider which uses its built-in relation model). This is important to be able to select features within a certain distance, which features intersects the point that was clicked on the map or features that intersect other features. At the current stage most intersection-tests are only done on bounding-box level, which - although very fast - isn't very precise.

NetTopologySuite is another great .NET library that among other does a lot of the things that SharpMap not (yet) can do. Fortunately Diego from the NTS team has made a link between NTS and SharpMap making this possible. They have posted an example on the NTS website showing how to check for exact intersection in SharpMap. Check it out. It is very cool and very powerful.

A basic survival guide to coordinate systems

Charlie Savage has a great post introducing coordinate systems in his blog. I too have had many questions on coordinate systems, projections and datums, and Charlie have made a great introduction to the world of "spatial reference systems". He promises to post some more on the topic, so keep track of his blog if you find these topics confusing (I know I did when I first was told about datums and coordinate systems).

Calculate the distance between two points on the earth

Here is a small C# method for calculating distances between two longitude/latitude points based on Movable Type's article on the subject. Since the earth isn't perfectly round, and since I've found several values for the approximate radius of the earth, this method is only approximate, but should be sufficiently accurate for most applications.

/// <summary>
/// Calculates the distance between to lat/long points and returns the approximate distance in kilometers
/// </summary>
/// <param name="from">Point in long/lat decimal degrees</param>
/// <param name="to">Point in long/lat decimal degrees</param>
/// <returns>Distance in kilometers</returns>
private double CalcDistance(Point from, Point to)
{
    double rad = 6371; //Earth radius in Km
    //Convert to radians
    double p1X = from.X / 180 * Math.PI;
    double p1Y = from.Y / 180 * Math.PI;
    double p2X = to.X / 180 * Math.PI;
    double p2Y = to.Y / 180 * Math.PI;
    
    return Math.Acos(Math.Sin(p1Y) * Math.Sin(p2Y) +
        Math.Cos(p1Y) * Math.Cos(p2Y) * Math.Cos(p2X - p1X)) * rad;
}

If you want the distance returned in something else than kilometers, change the radius to whatever unit you would like the output unit to be in.

I have successfully incorporated this in a SharpMap-based application for doing nearest-object searches in a map. It is fun to see how the earth-curvate is taken into account; -this is particularily visible towards the poles in a flat projected map.

Using SharpMap as a Google Earth Network Link

Since Google Earth is so popular, I thought why not join in and have SharpMap serve maps for GE? It is quite simple to make your own GE Network Link, and here are a few pointers, and you can download a demo web application you can use as well.

The GE Network Link basically works by requesting an XML webpage (well Google calls it KML), telling GE where to get the image tiles and where they should be placed on the globe. GE also adds a 'BBOX=minx,miny,maxx,maxy' querystring to the request so you can return a custom KML response containing the most appropriate tile(s). A response could look like this:

<?xml version="1.0" encoding="UTF-8" ?>
<kml xmlns="
http://earth.google.com/kml/2.0">
<Document>
  <open>1</open>
  <GroundOverlay>
    <open>1</open>
    <Icon>
      <href>http://localhost/TileServer.ashx?BBOX=-180,-90,0,90&size=512</href>
    </Icon>
    <LatLonBox>
      <north>90</north>
      <south>-90</south>
      <east>0</east>
      <west>-180</west>
      <rotation>0</rotation>
    </LatLonBox>
  </GroundOverlay>
</Document>
</kml>

You can add as many "GroundOverlays" as you like. For instance you can cover a requested BBOX with as many tiles as you like. That way you can increase the image size/resolution or cover a larger or smaller area than initially requested while getting quicker responses. You can read more about the KML format here.

Creating the KML is pretty straight-forward in ASP.NET, so I wont cover that here (but download the example and see for yourself).

Creating the image tile is easily done using SharpMap. Add an HttpHandler to you webpage, and create the following ProcessRequest method to render and return the map:

public void ProcessRequest(HttpContext context)

{

    //Get tilesize in request url

    int tilesize = int.Parse(context.Request.QueryString["size"]);

    //Get boundingbox requested

    string[] strbox = context.Request.QueryString["BBOX"].Split(new char[] { ',' });

    SharpMap.Geometries.BoundingBox box = new SharpMap.Geometries.BoundingBox

            (double.Parse(strbox[0]), double.Parse(strbox[1]),

            double.Parse(strbox[2]), double.Parse(strbox[3]));

 

    //Call custom method that sets up the map with the requested tilesize

    SharpMap.Map myMap = MapHelper.InitializeMap(new System.Drawing.Size(tilesize, tilesize));

    //Center on requested tile and set the appropriate view width

    myMap.Center = box.GetCentroid();

    myMap.Zoom = box.Width;

 

    //Render the map

    System.Drawing.Bitmap b = (System.Drawing.Bitmap)myMap.GetMap();

    //Create a PNG image which supports transparency

    context.Response.Clear();

    context.Response.Expires = 10000000;

    context.Response.ContentType = "image/png";

    System.IO.MemoryStream MS = new System.IO.MemoryStream();

    //Set background to the transparent color which will be see-through in Google Earth

    b.MakeTransparent(myMap.BackColor);

    b.Save(MS, System.Drawing.Imaging.ImageFormat.Png);

    byte[] buffer = MS.ToArray();

    //Send image response

    context.Response.OutputStream.Write(buffer, 0, buffer.Length);

    // tidy up 

    b.Dispose();

    context.Response.End();

}

The last thing we need to do is add the network link to GE. The easiest way to do this is to do this from within GE. From the menu Select Add –> Network link. Type a name for your network link and set the location to the URL of the XML service. Under “View-Based Refresh” , set the “When” parameter to “After camera stops”, and click OK, and you should be set to go !

Unfortunately GE doesn't provide you with the same smooth image transitions that it use for its own tiles, so you will have to do with the rather crude way of showing big red squares until the tiles have been loaded. Furthermore the BBOX GE requests isn't always very accurate, especially if you are looking south at an angle you will notice the BBOX is very much off.

The small demo application you can download below shows a world map with each country’s population density shown on a colorscale from blue to red. Copy the files to a virtual directory (make sure it runs as its own web-application, at least the sub-folders are located at the root of the webapplication). Set the network-link to point to the URL of the default.aspx page.

Download GoogleEarthServer.zip (270,65 KB)

GE also supports vectordata to be served as KML. It could be interesting to see who comes up with a "SharpMap KML Vector Server" first :-)

Adding SharpMap geometry to a PostGIS database

...or how to upload an ESRI Shapefile to PostGIS.

I've often been asked how to copy data from a shapefile to a PostGIS database. PostGIS comes with a commandline-tool for this (shp2pgsql.exe), but all it does is generate a huge SQL-textfile that you will need to run afterwards - Not very efficient I think - especially with the ASCII-representation of the geometry. Furthermore I've had several problems with it regarding many international letters in the attributes.

So why not try to let Npgsql and SharpMap do the job?

I've been working a bit with a small tool that makes it easy to upload an entire shapefile to a PostGreSQL/PostGIS database using SharpMap.

Below are some of the PostGIS/SharpMap related code explained:


 

 

First we create a Npgsql connection

NpgsqlConnection conn = new NpgsqlConnection("Server=localhost;Port=5432;User Id=username;Password=password;Database=myGisDB;")
NpgsqlCommand command = new NpgsqlCommand
();
command.Connection = conn;

The next step is to add a geometry column (see in the full source on how you create the table with all the attributes). In this case we set the spatial reference ID to '-1' and name the geometry column 'geom'.

command.CommandText = "SELECT AddGeometryColumn('','myTable','geom','-1','GEOMETRY',2);";
command.ExecuteNonQuery();

Now we are ready to upload to the database, so lets get hold of that shapefile! First we set up a datasource:

SharpMap.Data.Providers.ShapeFile shp = new SharpMap.Data.Providers.ShapeFile(@"C:\data\MyShape.shp", false);

We can now query all the feature object IDs, by using an extents-query on the full extents:

conn.Open();
List<uint> indexes = shp.GetObjectIDsInView(shp.GetExtents());

...and then loop through all the features:

foreach (uint idx in indexes)
{
   SharpMap.Data.FeatureDataRow
feature = shp.GetFeature(idx);
   
command.CommandText = "INSERT INTO \"myTable\" (\"geom\") VALUES (GeomFromWKB(:geom,:srid));"
;
   command.Parameters.Add(":geom", NpgsqlTypes.NpgsqlDbType
.Bytea);
   command.Parameters[":geom"
].Value = feature.Geometry.AsBinary(); //Add the geometry as Well-Known Binary
   command.Parameters.Add(":srid", NpgsqlTypes.NpgsqlDbType
.Integer);
   //Set the SRID of the geometry - this must match the SRID we used when the column was created
   command.Parameters[":srid"].Value = -1;

   
//TODO: Add parameters for remaining columns if nessesary (in that case alter the INSERT commandtext accordingly) 
 
   command.ExecuteNonQuery();
}
//Clean up
conn.Close();
shp.Close();

...and that is all there is to it !

The great thing about this, is that it is easy to change this to take any other SharpMap datasource and upload as well. And with Christians OGR extension you can suddenly upload a bunch of datasource directly to PostGIS.

Download the full source and compiled binaries here: Shape2Pgsql.zip (624,3 KB) (updated April 26, 2006)

Working with 8bit images in .NET

For some reason Microsoft made it pretty hard to work with anything but 24bit images. This is even though they provide you with several pixel formats, but no way of setting and getting the values of a pixel. If you use the SetPixel(x,y) or GetPixel(x,y) methods, your application will fail. There are plenty of articles and blogs on the Internet on how to do direct access on 1bit and 24bit, but I wasn't able to find anything on 8bit.

This article will cover some of the basics on how to access 8 bit greyscale or indexed images, by accessing the bitmapdata directly in memory. This also has the benefit of being much faster than the Set/GetPixel methods provided by the .NET Framework.

Before we can access the memory directly, we must lock its place in memory. We can do this by calling the Bitmap.LockBits() method:

BitmapData bmd  = myBitmap.LockBits(new Rectangle(0, 0, myBitmap.Width, myBitmap.Height),
ImageLockMode.ReadWrite, myBitmap.PixelFormat);

Likewise when we are done using the BitmapData, remember to unlock the data:

myBitmap.UnlockBits(bmd);

Now we need a method that can access the BitmapData. Lets make our own SetPixel and GetPixel method. Here we assume that we are dealing with 8bit pixels. We also add the 'unsafe' keyword since direct memory access isn't thread safe. I won't cover the Stride and Scan0 values. Bob Powell has a nice article on this.

public unsafe void SetPixel(int x, int y, byte c)
{
	byte* p = (byte *)bmd.Scan0.ToPointer();
	int offset=y*bmd.Stride+(x);
	p[offset] = c;
}

public unsafe Byte GetPixel(int x, int y)
{
	byte* p = (byte *)bmd.Scan0.ToPointer();
	int offset=y*bmd.Stride+x;
	return p[offset];
}

It is worth noting that GetPixel only returns a byte and not a color. The byte represents a number between 0 and 255. Each of the values is actually an index to a color palette. The palette could specify that for instance index 0 is black, index 1 is red, index 3 is blue etc. If you want a greyscale image, we can override the color palette. Let's set index 0 to black, index 255 to white, and linearly distribute the grayscale in between.

public static void SetGrayscalePalette(Bitmap b)
{
	ColorPalette pal = b.Palette;
	for(int i = 0; i < 256; i++)
	pal.Entries[i] = Color.FromArgb( 255, i, i, i );
	b.Palette = pal;
}

You can easily override this palette to specify other than grayscale images.

We can likewise create a function that can convert an index to a System.Drawing.Color. If you are working with a grayscale image, there is probably no need for this.

public System.Drawing.Color GetColorFromIndex(byte c)
{
	return = myBitmap.Palette.Entries[c];
}

Now let's put it all together into an easy-to-use 8bit image access class. Remember to allow unsafe code blocks before compiling.

using System;
using System.Drawing;
using System.Drawing.Imaging;
using System.Runtime.InteropServices;
namespace ImageProc
{
   /// <summary>
   /// Class used for direct memory access to 8bit grayscale images
   /// </summary>
   public class Image8Bit : IDisposable
   {
      private BitmapData bmd;
      private Bitmap b;
      /// <summary>
      /// Locks an 8bit image in memory for fast get/set pixel functions.
      /// Remember to Dispose object to release memory.
      /// </summary>
/// Bitmap reference public Image8Bit (Bitmap bitmap) { if(bitmap.PixelFormat!=System.Drawing.Imaging.PixelFormat.Format8bppIndexed) throw(new System.Exception("Invalid PixelFormat. 8 bit indexed required")); b = bitmap; //Store a private reference to the bitmap bmd = b.LockBits(new Rectangle(0, 0, b.Width, b.Height),
ImageLockMode.ReadWrite, b.PixelFormat); } /// <summary> /// Releases memory /// </summary> public void Dispose() { b.UnlockBits(bmd); } /// <summary> /// Gets color of an 8bit-pixel /// </summary> /// <param name="x">Row</param> /// <param name="y">Column</param> /// <returns>Color of pixel</returns> public unsafe System.Drawing.Color GetPixel(int x, int y) { byte* p = (byte *)bmd.Scan0.ToPointer(); //always assumes 8 bit per pixels int offset=y*bmd.Stride+x; return GetColorFromIndex(p[offset]); } /// <summary> /// Sets color of an 8bit-pixel /// </summary> /// <param name="x">Row</param> /// <param name="y">Column</param> /// <param name="c">Color index</param> public unsafe void SetPixel(int x, int y, byte c) { byte* p = (byte *)bmd.Scan0.ToPointer(); //always assumes 8 bit per pixels int offset=y*bmd.Stride+(x); p[offset] = c; } /// <summary> /// Sets the palette for the referenced image to Grayscale /// </summary> public void MakeGrayscale() { SetGrayscalePalette(this.b); } /// <summary> /// Sets the palette of an image to grayscales (0=black, 255=white) /// </summary> /// <param name="b">Bitmap to set palette on</param> public static void SetGrayscalePalette(Bitmap b) { ColorPalette pal = b.Palette; for(int i = 0; i < 256; i++) pal.Entries[i] = Color.FromArgb( 255, i, i, i ); b.Palette = pal; } private System.Drawing.Color GetColorFromIndex(byte c) { return = b.Palette.Entries[c]; } } }

Blog moved !

For a long time I wanted to start using a better blog-engine, and have now finally upgraded to dasBlog, a great ASP.NET-based blogging engine. At the same time I moved it to a new URL, so from now on, you can find my blog at www.sharpgis.net