shpreader.cs

实现SHP

///GeoCon, free tool to create gml & svg from gis files. 
///Copyright(C) 2005 Amri Rosyada
///Distributed under GNU-LGPL, see a copy of the license in root directory

using System;
using System.IO;
using System.Runtime.Serialization.Formatters.Binary;

using GeoCon;
using GeoCon.Data;
using GeoCon.dbf;
using gml;

namespace GeoCon.shp
{
	/// <summary>
	/// A class for reading arcview shapefile(*.shp), including its index(*.shx) and data(*.dbf).
	/// </summary>
	public class ShapeReader
	{
		#region private fields
		private short shpFileCode = 9994;
		private int RecordCount;
		private int currentShapeType;

		private BinaryReader shpReader;
		private BinaryReader shxReader;
		private dbf.DbfReader dbfReader;
		private string name;
		//private const string defaultsrs = "http://www.opengis.net/gml/srs/epsg.xml#4326";
		private string defaultsrs = " ";

		//private dbf.Header dbHeader;
		private shp.RecordHeader[] recordIndex;
		#endregion

		#region public fields

		/// <summary>
		/// Occurs when there's noticable change in this reader progress.
		/// </summary>
		public event GeoCon.Data.StatusEventHandler StatusChange;
		private GeoCon.Data.StatusEventArgs StatusArgs = new GeoCon.Data.StatusEventArgs("",0); 

		/// <summary>Holds various info about the shapefile</summary>
		public shp.ShapeHeader mainHeader;
		//public gml.AbstractGeometryType[] Shapes;

		/// <summary>Field collection of dbf and geometry data.</summary>
		public Data.FieldCollection Fields;
		#endregion

		#region constructor
		/// <summary>
		/// Constructs new ShapeReader
		/// </summary>
		/// <param name="shapepath">the full path to a shapefile(*.shp)</param>
		public ShapeReader(string shapepath)
		{
			try
			{
				string fne =System.IO.Path.GetFileNameWithoutExtension(shapepath);
				name=fne;
				string shxpath = fne + ".shx";
				string dbfpath = fne + ".dbf";

				shpReader = new BinaryReader(File.Open(shapepath,FileMode.Open,FileAccess.Read,FileShare.Read));
				shxReader = new BinaryReader(File.Open(shxpath,FileMode.Open,FileAccess.Read,FileShare.Read));
				if(File.Exists(dbfpath)) dbfReader = new dbf.DbfReader(dbfpath);
				shxpath=null;
				dbfpath=null;
			}
			catch(FileNotFoundException fnfe)
			{
				//System.Windows.Forms.MessageBox.Show("file '"+fnfe.FileName+"' not found.");
				throw new FileNotFoundException(fnfe.Message); 
			}
			catch(System.Exception ioe)
			{
				//System.Windows.Forms.MessageBox.Show("error : \r"+ioe.Message);
				throw new Exception(ioe.Message,ioe.InnerException);   
			}

		}
		#endregion

		#region public methods

		/// <summary>
		/// Begin read .shp, .shx, and .dbf file.
		/// </summary>
		/// <returns>The MapData which holds shapes & data.</returns>
		public MapData Read()
		{
			UpdateStatus("reading "+name+"..",0);
			MapData mapdata = new MapData();
			mapdata.Name=name;
			mainHeader = new shp.ShapeHeader();
			if(dbfReader!=null) dbfReader.Read(ref mapdata);
			readShapeHeader(ref mapdata);
			readIndexHeader();
			readRecordHeader(ref mapdata);
			mapdata.Info.FileDirectory = System.Environment.CurrentDirectory;//System.IO.Path.GetDirectoryName(name);
			mapdata.Info.FileGeometry  = name+".shp";
			mapdata.Info.FileData =(dbfReader!=null)? name+".dbf" : "not found";
			mapdata.Info.FileIndex = name+".shx";
			mapdata.Info.RecordsCount = RecordCount;

			UpdateStatus("Finished reading "+name,0);

			return mapdata;
		}
		/// <summary>
		/// close readers/streams used in this reader.
		/// </summary>
		public void Close()
		{	
			if(dbfReader!=null) dbfReader.Close();
			shpReader.Close();
			shxReader.Close();
		}
		#endregion

		#region read header
		/// <summary>
		/// 
		/// </summary>
		/// <param name="mapdata"></param>
		private void readShapeHeader(ref MapData mapdata)
		{
			shp.ShapeFileHeader fh = mainHeader.shape;
			mainHeader.shape.FileCode = SwapEndian32(shpReader.ReadInt32());
			if (mainHeader.shape.FileCode != shpFileCode)
			{
				shpReader.Close();
				throw new Exception("filetype unknown, "+mainHeader.shape.FileCode.ToString());
			}
			//skip 5 * int32, me know not what are they
			shpReader.BaseStream.Seek(20,SeekOrigin.Current);
			mainHeader.shape.FileLength = SwapEndian32(shpReader.ReadInt32());
			mainHeader.shape.version = shpReader.ReadInt32();
			mainHeader.shape.shapeType = shpReader.ReadInt32();

			mapdata.Info.FileVersion = mainHeader.shape.version.ToString();
			mapdata.Info.BoundsLeft = shpReader.ReadDouble();
			mapdata.Info.BoundsBottom = shpReader.ReadDouble();
			mapdata.Info.BoundsRight = shpReader.ReadDouble();
			mapdata.Info.BoundsTop = shpReader.ReadDouble();

			mapdata.BoundingBox  = new gml.BoxType(new double[4]{mapdata.Info.BoundsLeft,mapdata.Info.BoundsTop,mapdata.Info.BoundsRight,mapdata.Info.BoundsBottom});
		}

		/// <summary>
		/// 
		/// </summary>
		private void readIndexHeader()
		{
			//byte start offset = 1
			mainHeader.index.FileCode = SwapEndian32(shxReader.ReadInt32());

			//skip 5 of 4-bytes integer , me know not what them
			shxReader.BaseStream.Seek(20,SeekOrigin.Current);

			//byte start offset = 25
			mainHeader.index.FileLength = SwapEndian32(shxReader.ReadInt32());
			RecordCount = (mainHeader.index.FileLength - 50) / 4;
		}

		/// <summary>
		/// Reads header for each shape. 
		/// </summary>
		/// <param name="mapdata"></param>
		private void readRecordHeader(ref MapData mapdata)
		{
			recordIndex = new shp.RecordHeader[RecordCount];
			GeometryField fgeo = new Data.GeometryField();
			mapdata.Fields.Add(fgeo);
			
			shxReader.BaseStream.Seek(100,SeekOrigin.Begin);

			for (int i=0;i<recordIndex.Length;i++)
			{
				UpdateStatus("",100*i/RecordCount);

				recordIndex[i] = new shp.RecordHeader();
				recordIndex[i].indexOffset = SwapEndian32(shxReader.ReadInt32());
				recordIndex[i].indexLength = SwapEndian32(shxReader.ReadInt32());

				//TODO ?: check if shpoff bigger than the stream length ;
				long shpoff = recordIndex[i].indexOffset*2;	//we use long for large files;
				shpReader.BaseStream.Seek(shpoff,SeekOrigin.Begin);

				recordIndex[i].shapeNumber = SwapEndian32(shpReader.ReadInt32());
				recordIndex[i].shapeLength = SwapEndian32(shpReader.ReadInt32());

				//shape number start from 1
				if(recordIndex[i].shapeNumber!=(i+1)) throw new Exception("record number in shape file doesn't match index file : " + recordIndex[i].shapeNumber.ToString()+" vs "+(i+1).ToString());
				currentShapeType = shpReader.ReadInt32();
				if(currentShapeType != mainHeader.shape.shapeType) throw new Exception("shape type in record no "+i.ToString()+" doesn't match with the header");

				readSingleShape(i,ref fgeo);
			}
		}
		#endregion

		#region read shapes

		/// <summary>
		/// Read arcview shape
		/// </summary>
		/// <param name="index">record index, zero based</param>
		/// <param name="fgeo"></param>
		private void readSingleShape(int index, ref GeometryField fgeo)
		{
			if(index<recordIndex.GetLowerBound(0) || index > recordIndex.GetUpperBound(0) ) 
			{
				throw new ArgumentOutOfRangeException("index must be between "+recordIndex.GetLowerBound(0).ToString()+" and "+recordIndex.GetUpperBound(0).ToString()) ;
			}
			switch (currentShapeType)
			{
				case (int)shp.TYPE.typeNullShape :
					fgeo.Add(fgeo.NullSymbol);
					break;
				case (int)shp.TYPE.typePoint :
					readPoint(ref fgeo);
					break;
				case (int)shp.TYPE.typePointM :
					readPointM(ref fgeo);
					break;
				case (int)shp.TYPE.typePointZ :
					readPointM(ref fgeo);
					break;
				case (int)shp.TYPE.typeMultiPoint :
					readMultiPoint(ref fgeo);
					break;
				case (int)shp.TYPE.typeMultiPointM :
					readMultiPointM(ref fgeo);
					break;
				case (int)shp.TYPE.typeMultiPointZ :
					readMultiPointM(ref fgeo);
					break;
				case (int)shp.TYPE.typePolyLine :
					readPolyline(ref fgeo);
					break;
				case (int)shp.TYPE.typePolyLineM :
					readPolylineM(ref fgeo);
					break;
				case (int)shp.TYPE.typePolyLineZ :
					readPolylineM(ref fgeo);
					break;
				case (int)shp.TYPE.typePolygon :
					readPolygon(ref fgeo,false);
					break;
				case (int)shp.TYPE.typePolygonM :
					readPolygon(ref fgeo,true);
					break;
				case (int)shp.TYPE.typePolygonZ :
					readPolygon(ref fgeo,true);
					break;
				case (int)shp.TYPE.typeMultiPatch :
					readMultiPatch(ref fgeo,true);
					break;
				default:
					break;
			}
		}
		#endregion

		#region read point, pointM & Z
		/// <summary>
		/// Read a point
		/// </summary>
		/// <param name="fgeo"></param>
		private void readPoint(ref GeometryField fgeo)
		{
			double _x = shpReader.ReadDouble();
			double _y = shpReader.ReadDouble();
			gml.CoordinatesType _ct = new gml.CoordinatesType(new double[2]{_x,_y});
			gml.PointType pt = new gml.PointType(_ct);
			fgeo.Add(new gml.PointPropertyType(pt));
		}
		/// <summary>
		/// Reads PointM or PointZ, M values treated as Z values.
		/// </summary>
		/// <param name="fgeo"></param>
		private void readPointM(ref GeometryField fgeo)
		{
			double _x = shpReader.ReadDouble();
			double _y = shpReader.ReadDouble();
			double _z = shpReader.ReadDouble();
			gml.CoordinatesType _ct = new gml.CoordinatesType(new double[2]{_x,_y});
			_ct.hasZ = true;
			_ct.zArray=new double[1]{_z};
			gml.PointType pt = new gml.PointType(_ct);
			fgeo.Add(new gml.PointPropertyType(pt));
		}

		#endregion

		#region read multipoint, multipointM & Z
		/// <summary>
		/// Read multipoint
		/// </summary>
		/// <param name="fgeo"></param>
		private void readMultiPoint(ref GeometryField fgeo)
		{
			//multipoint bounding box
			double _left = shpReader.ReadDouble();
			double _bottom = shpReader.ReadDouble();
			double _right = shpReader.ReadDouble();
			double _top = shpReader.ReadDouble();

			int _partsCount = shpReader.ReadInt32();
			gml.PointType[] pts = new gml.PointType[_partsCount]; 
			double _x, _y;
			for (int j=0;j<_partsCount;j++)
			{
				_x = shpReader.ReadDouble();
				_y = shpReader.ReadDouble();
				gml.CoordinatesType _ct = new gml.CoordinatesType(new double[2]{_x,_y});
				pts[j] = new gml.PointType(_ct);
			}
			gml.MultiPointType mpt = new gml.MultiPointType(defaultsrs,pts);
			fgeo.Add(new gml.MultiPointPropertyType(mpt)); 
			pts=null;
			mpt=null;
		}
		/// <summary>
		/// Read MultiPointM or MultiPointZ, M values treated as Z values
		/// </summary>
		/// <param name="fgeo"></param>
		private void readMultiPointM(ref GeometryField fgeo)
		{
			//multipoint bounding box
			double _left = shpReader.ReadDouble();
			double _bottom = shpReader.ReadDouble();
			double _right = shpReader.ReadDouble();
			double _top = shpReader.ReadDouble();

			int _partsCount = shpReader.ReadInt32();
			gml.CoordinatesType[] cts=new CoordinatesType[_partsCount];
			double _x, _y;
			for (int j=0;j<_partsCount;j++)
			{
				_x = shpReader.ReadDouble();
				_y = shpReader.ReadDouble();
				cts[j]=new gml.CoordinatesType(new double[2]{_x,_y});
			}
			double miniValue = shpReader.ReadDouble();
			double maxiValue = shpReader.ReadDouble();
			gml.PointType[] pts = new gml.PointType[_partsCount]; 
			for (int j=0;j<_partsCount;j++) //each part contains one point
			{
				cts[j].hasZ=true;
				cts[j].zArray=new double[1]{shpReader.ReadDouble()};
				pts[j]=new PointType(cts[j]);
			}
			if(_partsCount>1)
			{
				gml.MultiPointType mpt = new gml.MultiPointType(defaultsrs,pts);
				fgeo.Add(new gml.MultiPointPropertyType(mpt)); 
				mpt=null;
			}
			else
			{
				fgeo.Add(new gml.PointPropertyType(pts[0])); 
			}
			cts=null;
			pts=null;
		}

		#endregion

		#region read polyline, polylineM & Z
		/// <summary>
		/// Reads a polyline
		/// </summary>
		/// <param name="fgeo"></param>
		private void readPolyline(ref GeometryField fgeo)
		{
			double _left = shpReader.ReadDouble();
			double _bottom = shpReader.ReadDouble();
			double _right = shpReader.ReadDouble();
			double _top = shpReader.ReadDouble();

			int _partsCount = shpReader.ReadInt32();
			int _totalPoints = shpReader.ReadInt32();
			int kend, dcount;
			gml.LineStringType[] LSArray = new gml.LineStringType[_partsCount]; 
			int[] _partPoints = readIntegerArray(shpReader,_partsCount);
			for (int j=0;j<_partPoints.Length;j++)
			{
				kend = (j!=(_partsCount-1))? _partPoints[j+1] : _totalPoints;
				dcount = (kend-_partPoints[j]);
				gml.CoordinatesType _ct = new gml.CoordinatesType(ReadTuples(shpReader,dcount));
				LSArray[j] = new gml.LineStringType(_ct); 
			}
			if(_partsCount>1) //multilinestring
			{
				gml.MultiLineStringType mls = new gml.MultiLineStringType(defaultsrs,LSArray);
				fgeo.Add(new MultiLineStringPropertyType(mls));
			}
			else //linestring
			{
				fgeo.Add(new gml.LineStringPropertyType(LSArray[0]));
			}
			LSArray=null;
			_partPoints=null;
		}
		/// <summary>
		/// Reads a polylineM or polylineZ, M values treated as Z values.
		/// </summary>
		/// <param name="fgeo">GeometryField to store the shape</param>
		private void readPolylineM(ref GeometryField fgeo)
		{
			double _left = shpReader.ReadDouble();
			double _bottom = shpReader.ReadDouble();
			double _right = shpReader.ReadDouble();
			double _top = shpReader.ReadDouble();

			int _partsCount = shpReader.ReadInt32();
			int _totalPoints = shpReader.ReadInt32();
			int kend, dcount;
			gml.CoordinatesType[] cts = new CoordinatesType[_partsCount];
			int[] _partPoints = readIntegerArray(shpReader,_partsCount);
			for (int j=0;j<_partPoints.Length;j++)
			{
				kend = (j!=(_partsCount-1))? _partPoints[j+1] : _totalPoints;
				dcount = (kend-_partPoints[j]);
				cts[j]=new gml.CoordinatesType(ReadTuples(shpReader,dcount));
			}
			double miniValue = shpReader.ReadDouble(); //or startlength?
			double maxiValue = shpReader.ReadDouble(); //or endlength?
			for (int j=0;j<_partPoints.Length;j++)
			{
				kend = (j!=(_partsCount-1))? _partPoints[j+1] : _totalPoints;
				dcount = (kend-_partPoints[j]);
				cts[j].hasZ=true;