Using the NetTopologySuite to read and write Shapefiles in C#

I was looking for a .NET library to use to read and write Shapefiles without using unmanaged dynamic link libraries. I found quite a few that could read but not write Shapefiles. In the end, I found the NetTopologySuite at http://code.google.com/p/nettopologysuite/ to suit my requirements the best.

The following code snippet show how to use the library's ShapefileDataReader to read a Shapefile. In the example, all the database records and the Shape geometries are read into a features collection.

using GisSharpBlog.NetTopologySuite.Features;

using GisSharpBlog.NetTopologySuite.Geometries;

using GisSharpBlog.NetTopologySuite.IO;

using GeoAPI.Geometries;

...etc....

GeometryFactory factory = new GeometryFactory();

ShapefileDataReader shapeFileDataReader = new ShapefileDataReader(shpFilename, factory);

//Display the shapefile type

ShapefileHeader shpHeader = shapeFileDataReader.ShapeHeader;

Console.WriteLine(string.Format("Shape type: {0}", shpHeader.ShapeType));

//Display the min and max bounds of the shapefile

IEnvelope bounds = shpHeader.Bounds;

Console.WriteLine(string.Format("Min bounds: ({0},{1})", bounds.MinX, bounds.MinY));

Console.WriteLine(string.Format("Max bounds: ({0},{1})", bounds.MaxX, bounds.MaxY));

//Display summary information about the Dbase file

DbaseFileHeader header = shapeFileDataReader.DbaseHeader;

Console.WriteLine("Dbase info");

Console.WriteLine(string.Format("{0} Columns, {1} Records", header.Fields.Length, header.NumRecords));

for (int i = 0; i < header.NumFields; i++)

{

DbaseFieldDescriptor fldDescriptor = header.Fields[i];

Console.WriteLine(string.Format("   {0} {1}", fldDescriptor.Name, fldDescriptor.DbaseType));

}

//Reset the pointer to the start of the shapefile, just in case

//shapeFileDataReader.Reset();

//Read through all records of the shapefile (geometry and attributes) into a feature collection 

ArrayList features = new ArrayList();

while (shapeFileDataReader.Read())

{

Feature feature = new Feature();

AttributesTable attributesTable = new AttributesTable();

string[] keys = new string[header.NumFields];

IGeometry geometry = (Geometry)shapeFileDataReader.Geometry;

for (int i = 0; i < header.NumFields; i++)

{

DbaseFieldDescriptor fldDescriptor = header.Fields[i];

keys[i] = fldDescriptor.Name;

attributesTable.AddAttribute(fldDescriptor.Name, shapeFileDataReader.GetValue(i));

}

feature.Geometry = geometry;

feature.Attributes = attributesTable;

features.Add(feature);

}

//Close and free up any resources

shapeFileDataReader.Close();

shapeFileDataReader.Dispose();

It is important to free up any resources used by calling the Close and Dispose methods, especially if you want to apply any file based operations e.g. deleting or renaming onto the Shapefile.

The following screenshot shows an example output from the code snippet.


The ShapefileDataWriter class can be used to create a Shapefile. It seems that the class can only create new Shapefiles and there are no methods append to an existing Shapefile. The following code snippet shows how to write features with geometry and attributes into a Shapefile.

//Create a new shapefile from features

GeometryFactory outGeomFactory = new GeometryFactory();

//The constructor will append a ".shp" extension to the input filename so the 1st argument should not have an extension at all

string outShpFilenameNoExt = Path.Combine(Path.GetDirectoryName(outShpFilename), Path.GetFileNameWithoutExtension(outShpFilename));

ShapefileDataWriter writer = new ShapefileDataWriter(outShpFilenameNoExt, outGeomFactory);

DbaseFileHeader outDbaseHeader = ShapefileDataWriter.GetHeader((Feature)features[0], features.Count);

writer.Header = outDbaseHeader;

writer.Write(features); 

WebApp to show DBASE (*.dbf) file information

The ESRI Shapefile format stores database attributes in DBASE DBF (*.dbf) file and shape geometries in another associated file (*.shp). Using the DBASE file structure documentation in http://www.dbf2002.com/dbf-file-format.html, I wrote a simple Javascript web app to show basic information about a DBF file, similar to the open source ShapeLib's dbfinfo executable.

 The web app can be run from this page http://dominoc925-pages.appspot.com/webapp/dbf_info/default.html.

1. To use the web app, simply drag and drop one or more DBASE files (*.dbf) into the dashed box as shown in the screenshot below.
   
   
2. Or click the button and choose one or more (*.dbf) files. .
   
   
   
   Basic information about the DBF file is displayed, including the number of records, the number of columns, column names and type.
   
   
   
   

Note: This will run only on Chrome and a development version of FireFox that uses Gecko 7 at the moment. 

WebApp to show Shapefile information

I wrote a simple HTML5 Web App to show basic information about one or more ESRI Shapefile (*.shp) files similar to the shpinfo command from the open source ShapeLib library. The processing of the Shapefile is done in the local web browser without having to upload it to a server. The web app can be run from this web site http://dominoc925-pages.appspot.com/webapp/shpfile_info/default.html.

One the page is displayed, either click the button and select one or more Shapefiles or drag and drop *.shp files into the dashed box as shown in the screenshots below.

The Shapefile information is displayed in the web page.

Applying coordinate transformations to ESRI world files

The ESRI World file is a plain text file with numbers on six lines. These numbers actually are just the coefficients of an Affine transformation matrix. The world file helps to place the raster image on a 2-dimensional plane. Since the world file is just a matrix, it is a simple matter to perform coordinate transformations (translation, rotation, and scaling) to the raster image by applying the various transformation matrices with the world file matrix.

An example of a world file is shown below.

0.124966541755889
-0.216448399567377
-0.216438882466538
-0.124961047055157
2132.25596524947
1263.76014101498

If we label the lines as the following for discussion:
A
D
B
E
C
F

Then the Affine transformation matrix of the raster image is

A C E
B D F
0 0 1

For more information on the world file format, visit the Wikipedia site http://en.wikipedia.org/wiki/World_file.

Translation transformation
To apply a translation transformation to the raster image, it is just a matrix multiplication of the translation matrix with the world file matrix, as shown below. dx and dy are the x offset and y offset of the translation.

| 1 0 dx |   | A C E |
| 0 1 dy | * | B D F |
| 0 0 1  |   | 0 0 1 |

Scaling transformation
To apply a scaling transformation to the raster image, just multiply the world file matrix with the scaling transformation matrix. s is the scale factor.

Note: if the scaling origin is not at (0,0), then apply the translation transformation to move the origin to the scaling origin (x,y) first. After the scaling is done, the apply an inverse translation translation to move the origin back.

| s 0 0 |   | A C E |
| 0 s 0 | * | B D F |
| 0 0 1 |   | 0 0 1 |

Rotation transformation
Similarly, to apply a rotation transformation to the raster image, multiply the rotation transformation matrix with the world file matrix. r is the rotation angle in radians (positive clockwise). 

Note: if the rotation origin is not at (0,0), then apply the translation transformation to move the origin to the rotation origin (x,y) first. After the rotation is done, the apply an inverse translation translation to move the origin back.

|  cos(r) sin(r) 0 |   | A C E |
| -sin(r) cos(r) 0 | * | B D F |
|  0      0      1 |   | 0 0 1 |

C# code to determine if the polygon vertices are in the clockwise direction

In some GIS polygon data formats such as the ESRI shape file format, the holes of a polygon are stored in the counter-clockwise direction while the outer polygon is stored in the clockwise direction. If you are writing a program to read a polygon geometry containing holes, it may be necessary to determine the holes by checking whether the vertices are in the counter-clockwise or clockwise direction.

The C# code snippet below is useful to determine whether the vertices of a polygon are in the clockwise direction or counter clockwise direction. Simply pass in the polygon vertices into the function as an array of PointF structures where the first and last member of the array are the same point. The function will return true if the vertices are in the clockwise direction and false if they are in the counter-clockwise direction.

private bool IsClockwisePolygon(PointF[] polygon)

{

bool isClockwise = false;

double sum = 0;

for ( int i = 0; i < polygon.Length-1; i++)

{

sum += (polygon[i + 1].X - polygon[i].X) * (polygon[i + 1].Y + polygon[i].Y);

}

isClockwise = (sum > 0) ? true : false;

return isClockwise;

}

Google Gadget for creating ESRI World Files

Note: A WebApp version of this Google Gadget will be up in a short while to replace this gadget.

The ESRI World file is a text file used to georeference a corresponding image file. A detailed description of the format can be found in this Wikipedia page http://en.wikipedia.org/wiki/World_file.

Recently I had to create ESRI World files for hundreds of aerial photographs, each having different origin, resolution and orientation. As a result of that exercise, I learnt enough about the World File format to create this Google Gadget.

This gadget can help to create the world files given the following input parameters:

1. X coordinate of the upper left pixel in ground units,
2. Y coordinate of the upper left pixel in ground units,
3. Size of the pixel in the X-axis direction in ground units,
4. Size of the pixel in the Y-axis direction in ground units,
5. Rotation angle in degrees of the image from North (clockwise is positive) 

Just fill in the required input, click the button and copy and paste the results into a world file, which must be named according to the type of image file. For example, the World file should have the file extension of .tfw for associating with TIFF images.