Shapefile Z values "Elevation"

Topics: SharpMap Project, SharpMap v2.0
Aug 21, 2014 at 2:54 PM
Edited Aug 21, 2014 at 2:54 PM
Hi all,

I'm trying to extract the Z data from some Linestrings in a shapefile. So far I can successfully iterate through the attribute tables, and extract data per segment, but when I pull out my vertex coordinates for c.X,c.Y,c.Z the Z values give a non-number NaN.

From what I've read around and regarding the shapefile specification, the Z values are sometimes stored elsewhere in the file. What container would I then look for? I created the data myself using GlobalMapper v14.
Coordinator
Aug 21, 2014 at 6:45 PM
SharpMap's Shapefile provider disregards z- and m-ordinate values.
Aug 22, 2014 at 8:31 AM
Edited Aug 22, 2014 at 8:36 AM
Hi,
I also needed the 3d values for the shape files. Using the Esri whitepaper it was not too difficult to add the Z-values.
geoAPI.Geometries.Coordinate already has a z-value, but it is never used in sharpmap. If you succeed in filling it in, it stays available.
Be aware that all functionality within sharpmap disregards this z-value (for example the envelope will not take z into account).
Having said this, I only needed to change the ParseGeometry function to read the z-values for the shapefile and add them to the Coordinate objects.
I used this for the last month or so, but it is not yet tested on more complex elements. It is only tested on simple points (not tested on MultiPoint), simple polylines (each polyline just one part) and simple polygons (one outer ring, no inner rings). If someone can provide me with shape files containing more complex elements, I could test it on that.
Because there is nothing available to store a m-value, and because I did not need it, I left that alone. Anyone uses these ?

Here is my code for ParseGeometry
        private static IGeometry ParseGeometry(IGeometryFactory factory, ShapeType shapeType, BinaryReader brGeometryStream)
        {
            //Skip record number and content length
            brGeometryStream.BaseStream.Seek(8, SeekOrigin.Current); 

            var type = (ShapeType)brGeometryStream.ReadInt32(); //Shape type
                if (type == ShapeType.Null)
                    return null;

                if (shapeType == ShapeType.Point || shapeType == ShapeType.PointM )
                {
                    return factory.CreatePoint(new Coordinate(brGeometryStream.ReadDouble(), brGeometryStream.ReadDouble(),0));
                }
                if (shapeType == ShapeType.PointZ)
                {
                    return factory.CreatePoint(new Coordinate(brGeometryStream.ReadDouble(), brGeometryStream.ReadDouble(), brGeometryStream.ReadDouble()));
                }

                if (shapeType == ShapeType.Multipoint || shapeType == ShapeType.MultiPointM ||
                    shapeType == ShapeType.MultiPointZ)
                {
                    brGeometryStream.BaseStream.Seek(32, SeekOrigin.Current); //skip min/max box
                    var nPoints = brGeometryStream.ReadInt32(); // get the number of points
                    if (nPoints == 0)
                        return null;
                    var feature = new Coordinate[nPoints];
                    for (var i = 0; i < nPoints; i++)
                        feature[i] = new Coordinate(brGeometryStream.ReadDouble(), brGeometryStream.ReadDouble(),0);
                    if (shapeType == ShapeType.MultiPointZ)
                    {
                        brGeometryStream.BaseStream.Seek(16, SeekOrigin.Current); //skip min/max for Z
                        for (var i = 0; i < nPoints; i++)
                            feature[i].Z = brGeometryStream.ReadDouble();
                    }
                    return factory.CreateMultiPoint(feature);
                }

                if (shapeType == ShapeType.PolyLine || shapeType == ShapeType.Polygon ||
                    shapeType == ShapeType.PolyLineM || shapeType == ShapeType.PolygonM ||
                    shapeType == ShapeType.PolyLineZ || shapeType == ShapeType.PolygonZ)
                {
                    brGeometryStream.BaseStream.Seek(32, SeekOrigin.Current); //skip min/max box

                    var nParts = brGeometryStream.ReadInt32(); // get number of parts (segments)
                    if (nParts == 0 || nParts < 0)
                        return null;

                    var nPoints = brGeometryStream.ReadInt32(); // get number of points
                    var segments = new int[nParts + 1];
                    //Read in the segment indexes
                    for (var b = 0; b < nParts; b++)
                        segments[b] = brGeometryStream.ReadInt32();
                    //add end point
                    segments[nParts] = nPoints;

                    if ((int)shapeType % 10 == 3)
                    {
                        var lineStrings = new ILineString[nParts];
                        for (var lineID = 0; lineID < nParts; lineID++)
                        {
                            var line = new Coordinate[segments[lineID + 1] - segments[lineID]];
                            var offset = segments[lineID];
                            for (var i = segments[lineID]; i < segments[lineID + 1]; i++)
                                line[i - offset] = new Coordinate(brGeometryStream.ReadDouble(), brGeometryStream.ReadDouble(),0);
                            lineStrings[lineID] = factory.CreateLineString(line);
                        }
                        if (shapeType == ShapeType.PolyLineZ)
                        {
                            brGeometryStream.BaseStream.Seek(16, SeekOrigin.Current); //skip min/max for Z

                            for (var lineID = 0; lineID < nParts; lineID++)
                            {
                                var offset = segments[lineID];
                                for (var i = segments[lineID]; i < segments[lineID + 1]; i++)
                                    lineStrings[lineID].Coordinates[i - offset].Z = brGeometryStream.ReadDouble();
                            }
                        }

                        if (lineStrings.Length == 1)
                            return lineStrings[0];

                        return factory.CreateMultiLineString(lineStrings);
                    }

                    //First read all the rings
                    var rings = new ILinearRing[nParts];
                    for (var ringID = 0; ringID < nParts; ringID++)
                    {
                        var ring = new Coordinate[segments[ringID + 1] - segments[ringID]];
                        var offset = segments[ringID];
                        for (var i = segments[ringID]; i < segments[ringID + 1]; i++)
                            ring[i - offset] = new Coordinate(brGeometryStream.ReadDouble(), brGeometryStream.ReadDouble(),0);
                        rings[ringID] = factory.CreateLinearRing(ring);
                    }
                    if (shapeType == ShapeType.PolygonZ)
                    {
                        brGeometryStream.BaseStream.Seek(16, SeekOrigin.Current); //skip min/max for Z
                        for (var ringID = 0; ringID < nParts; ringID++)
                        {
                        var offset = segments[ringID];
                        for (var i = segments[ringID]; i < segments[ringID + 1]; i++)
                            rings[ringID].Coordinates[i - offset].Z = brGeometryStream.ReadDouble();
                        }
                    }

                    ILinearRing exteriorRing;
                    var isCounterClockWise = new bool[rings.Length];
                    var polygonCount = 0;
                    for (var i = 0; i < rings.Length; i++)
                    {
                        isCounterClockWise[i] = rings[i].IsCCW();
                        if (!isCounterClockWise[i])
                            polygonCount++;
                    }
                    if (polygonCount == 1) //We only have one polygon
                    {
                        exteriorRing = rings[0];
                        ILinearRing[] interiorRings = null;
                        if (rings.Length > 1)
                        {
                            interiorRings = new ILinearRing[rings.Length - 1];
                            Array.Copy(rings, 1, interiorRings, 0, interiorRings.Length);
                        }
                        return factory.CreatePolygon(exteriorRing, interiorRings);
                    }

                    var polygons = new List<IPolygon>();
                    exteriorRing = rings[0];
                    var holes = new List<ILinearRing>();

                    for (var i = 1; i < rings.Length; i++)
                    {
                        if (!isCounterClockWise[i])
                        {
                            polygons.Add(factory.CreatePolygon(exteriorRing, holes.ToArray()));
                            holes.Clear();
                            exteriorRing = rings[i];
                        }
                        else
                            holes.Add(rings[i]);
                    }
                    polygons.Add(factory.CreatePolygon(exteriorRing, holes.ToArray()));

                    return factory.CreateMultiPolygon(polygons.ToArray());
                }

            throw (new ApplicationException("Shapefile type " + shapeType + " not supported"));
        }
Aug 22, 2014 at 8:48 AM
A small extra note on the code:
It is mainly the same code as the 2D code, the only addition is an extra loop going over the z-values. The loop is inside an if statement (for example if (shapeType == ShapeType.PolyLineZ) )
This loop sets the z value of the coordinate elements already created in the XY loop. In the shapefile format, for each object you first have the XY,XY,XY... values for each vertex, then, if available in the same order the Z,Z,..., and finally M,M,M...
So that makes it easy, the loop going over the Z values is the same as the initial one going over the XY values.
Coordinator
Aug 22, 2014 at 9:13 AM
Nice effort. If you need the m-values also, you need to use some ICoordinateSequence implementation that can handle m-ordinates.
But that would be a different effort.

You can use NTSProvider and NTS Shapefile reading writing capabilities.
Aug 22, 2014 at 10:44 AM
FObermaier wrote:
You can use NTSProvider and NTS Shapefile reading writing capabilities.
Can you explain this ? I just changed the ParseGeometry function from the ShapeFile provider. So what would be the use of the NTSPovider ? Or is this remark related to the m-values ?
Coordinator
Aug 25, 2014 at 8:16 AM
I was referring to a case where you'd want to have the m-values, too.