shpreader.cs

实现SHP

				cts[j].zArray=ReadDoubles(shpReader,dcount);
			}
			if(_partsCount>1) //multilinestring
			{
				gml.LineStringType[] lss = new LineStringType[_partsCount];
				for(int i=0;i<cts.Length;i++) lss[i]=new LineStringType(cts[i]);
				gml.MultiLineStringType mls = new gml.MultiLineStringType(defaultsrs,lss);
				fgeo.Add(new MultiLineStringPropertyType(mls));
			}
			else //linestring
			{
				fgeo.Add(new gml.LineStringPropertyType(new gml.LineStringType(cts[0])));
			}
			_partPoints=null;
		}
		#endregion

		#region read polygon
	
		/// <summary>
		/// read polygon / multipolygon and add it into geometry field
		/// For better result: in Arcview, explode the multipart polygon first and then union the resulting polygons back ; plus do some cleaning.
		/// (in avenue using Explode,ReturnUnion,Clean) Must check if such script exists or write one.
		/// </summary>
		/// <param name="fgeo">geometry field where the polygon property will be added to.</param>
		/// <param name="hasMeasure">whether we are reading plain Polygon or PolygonM/PolygonZ.</param>
		private void readPolygon(ref GeometryField fgeo, bool hasMeasure)
		{
			//polygon bounding box
			double _left = shpReader.ReadDouble();
			double _bottom = shpReader.ReadDouble();
			double _right = shpReader.ReadDouble();
			double _top = shpReader.ReadDouble();

			int _partsCount = shpReader.ReadInt32();
			if(_partsCount<=0)
			{
				throw new System.Exception("encounter polygon with zero parts count"); 
			}
			else if(_partsCount==1) //singlepartpolygon
			{
				ReadSinglePartPolygon(ref fgeo,hasMeasure);
			}
			else //multipartpolygon
			{
				ReadMultiPartPolygon(ref fgeo,_partsCount,hasMeasure);
			}
		}

		/// <summary>
		/// Reads singlepart polygon
		/// </summary>
		/// <param name="fgeo">geometry field in which the polygonproperty will be added to</param>
		/// <param name="hasMeasure">whether we are reading plain Polygon or PolygonM/PolygonZ.</param>
		private void ReadSinglePartPolygon(ref GeometryField fgeo, bool hasMeasure)
		{
			int _totalPoints = shpReader.ReadInt32();
			int _partPoints = shpReader.ReadInt32(); //= _totalPoints, but we must advance anyway
			gml.CoordinatesType ct = new gml.CoordinatesType(ReadTuples(shpReader,_totalPoints));
			if(hasMeasure)
			{
				double miniValue = shpReader.ReadDouble();
				double maxiValue = shpReader.ReadDouble();

				ct.hasZ=true;
				ct.zArray=ReadDoubles(shpReader,_totalPoints);
			}
			gml.PolygonType pt = new gml.PolygonType(ct);
			fgeo.Add(new gml.PolygonPropertyType(pt));
		}

		/// <summary>
		/// Reads polygon with 2 or more parts.
		/// </summary>
		/// <param name="fgeo">geometry field in which the polygonproperty or multipolygonproperty will be added to</param>
		/// <param name="partsCount">number of parts</param>
		/// <param name="hasMeasure">whether we are reading plain Polygon or PolygonM/PolygonZ.</param>
		private void ReadMultiPartPolygon(ref GeometryField fgeo,int partsCount, bool hasMeasure)
		{
			int _totalPoints = shpReader.ReadInt32();
			int[] _partPoints = readIntegerArray(shpReader,partsCount);
			gml.BoxType[] boxes = new BoxType[partsCount];
			gml.CoordinatesType[] cts = new CoordinatesType[partsCount];
			for (int j=0;j<partsCount;j++)
			{	
				int dcount = (j != (partsCount-1))? (_partPoints[j+1] - _partPoints[j]) : (_totalPoints - _partPoints[j]);
				boxes[j]=new BoxType();
				cts[j] = new CoordinatesType( ReadTuplesWithBox(shpReader,dcount,ref boxes[j]) );
			}

			if(hasMeasure)
			{
				double miniValue = shpReader.ReadDouble();
				double maxiValue = shpReader.ReadDouble();
				for (int j=0;j<partsCount;j++)
				{	
					int dcount = (j != (partsCount-1))? (_partPoints[j+1] - _partPoints[j]) : (_totalPoints - _partPoints[j]);
					cts[j].hasZ=true;
					cts[j].zArray = ReadDoubles(shpReader,dcount);
				}
			}

			mbox[] ms = matchBoxes(boxes);
			int polycount=0;
			for(int i=0;i<ms.Length;i++)
			{
				if(ms[i].container!=-1) continue;
				polycount++;
				for(int j=0;j<ms.Length;j++) if(ms[j].container==i) ms[i].innerCount++;
			}
			gml.PolygonType[] pgs = new PolygonType[polycount];

			int polyindex=-1;
			for(int i=0;i<ms.Length;i++)
			{
				if(ms[i].container!=-1) continue;
				polyindex++;
				pgs[polyindex] = new PolygonType(cts[ms[i].index]);
				int icount=ms[i].innerCount;
				gml.PolygonTypeInnerBoundaryIs[] ibs = new PolygonTypeInnerBoundaryIs[icount];
				int innerindex=-1;
				for(int j=0;j<ms.Length;j++)
				{
					if(ms[j].container!=i) continue;
					innerindex++;
					ibs[innerindex]=new PolygonTypeInnerBoundaryIs(cts[ms[j].index]);
				}
				pgs[polyindex].innerBoundaryIs = ibs;
				ibs=null;
			}

			if(polycount==1)
			{
				fgeo.Add(new gml.PolygonPropertyType(pgs[0]));
			}
			else 
			{
				fgeo.Add(new gml.MultiPolygonPropertyType(new MultiPolygonType(defaultsrs,pgs)));
			}

			ms  = null;
			pgs = null;
			_partPoints = null;
			boxes = null;
		}

		/// <summary>
		/// Find relationships among parts in shape
		/// </summary>
		/// <param name="boxes">bounding box of each part </param>
		/// <returns>the containment relationship</returns>
		private mbox[] matchBoxes(gml.BoxType[] boxes)
		{
			mbox[] ms = new mbox[boxes.Length];
			for(int i=0;i<boxes.Length;i++)
			{ 
				ms[i]=new mbox();
				ms[i].index=i; 
				ms[i].container=-1;
				ms[i].innerCount=0;
			}
			for(int i=0;i<boxes.Length;i++)
			{
				for(int j=0;j<boxes.Length;j++)
				{
					if(i==j) continue;
					if(boxes[i].isContainedIn(boxes[j]))
					{
						if(ms[i].container==-1)
						{
							ms[i].container=j;
						}
						else
						{
							if(j==ms[i].container) continue;
							if(boxes[j].isContainedIn(boxes[ms[i].container]))	//nested parts
							{
								ms[j].container = ms[i].container; //immediate container
								ms[i].container = -1; //set the container's container to nothing --> new polygon;
								//TODO : Check if this algorithm gives correct result.., so far it's ok.
							}
						}
					}
				}
			}
			return ms;
		}

		/// <summary>
		/// private struct for storing relationship between parts
		/// </summary>
		private struct mbox
		{
			public int index;		//this shape index
			public int container;	//the container index
			public int innerCount;  //number of innerboundary
		}
		// Notes on poygon : 
		// 1. how do we determine if a part is outer or inner ?
		//    they're not always stored in fixed order
		//    they're not always 'clean'(no self intersection & it's inner boundary - there's if any - traverse counter-clockwise)
		// 2. how do we deal with inner boundary containing another inner boundary (island shapes)? 
		//    gml 2.0 innerboundary may not contain another boundary
		#endregion

		#region read multipatch
		/// <summary>
		/// read multipatch
		/// </summary>
		/// <param name="fgeo">Geometry field</param>
		/// <param name="hasMeasure">it must be true</param>
		private void readMultiPatch(ref GeometryField fgeo, bool hasMeasure)
		{
			//multipatch bounding box
			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[] _partPoints = readIntegerArray(shpReader,_partsCount);
			int[] _partTypes = readIntegerArray(shpReader,_partsCount);

			gml.BoxType[] boxes = new BoxType[_partsCount];
			gml.CoordinatesType[] cts = new CoordinatesType[_partsCount];
			for (int j=0;j<_partsCount;j++)
			{	
				int dcount = (j != (_partsCount-1))? (_partPoints[j+1] - _partPoints[j]) : (_totalPoints - _partPoints[j]);
				boxes[j]=new BoxType();
				if(_partTypes[j]==(int)PATCHTYPE.TRISTRIP || _partTypes[j]==(int)PATCHTYPE.TRIFAN)
				{
					double[] dd = ConvertToTriangles( ReadTuplesWithBox(shpReader,dcount,ref boxes[j]),_partTypes[j] );
					cts[j] = new CoordinatesType(dd);
				}
				else
				{
					cts[j] = new CoordinatesType( ReadTuplesWithBox(shpReader,dcount,ref boxes[j]) );
				}
			}
			if(hasMeasure)
			{
				double miniValue = shpReader.ReadDouble();
				double maxiValue = shpReader.ReadDouble();
				for (int j=0;j<_partsCount;j++)
				{	
					int dcount = (j != (_partsCount-1))? (_partPoints[j+1] - _partPoints[j]) : (_totalPoints - _partPoints[j]);
					cts[j].hasZ=true;
					if(_partTypes[j]==(int)PATCHTYPE.TRISTRIP || _partTypes[j]==(int)PATCHTYPE.TRIFAN)
					{
						cts[j].zArray = ConvertToTriangles( ReadDoubles(shpReader,dcount),_partTypes[j] );
					}
					else
					{
						cts[j].zArray = ReadDoubles(shpReader,dcount);
					}
					
				}
			}

			mbox[] ms = matchBoxes(boxes);
			int polycount=0;
			for(int i=0;i<ms.Length;i++)
			{
				if(ms[i].container!=-1) continue;
				polycount++;
				for(int j=0;j<ms.Length;j++) if(ms[j].container==i) ms[i].innerCount++;
			}
			gml.PolygonType[] pgs = new PolygonType[polycount];

			int polyindex=-1;
			for(int i=0;i<ms.Length;i++)
			{
				if(ms[i].container!=-1) continue;
				polyindex++;
				pgs[polyindex] = new PolygonType(cts[ms[i].index]);
				int icount=ms[i].innerCount;
				gml.PolygonTypeInnerBoundaryIs[] ibs = new PolygonTypeInnerBoundaryIs[icount];
				int innerindex=-1;
				for(int j=0;j<ms.Length;j++)
				{
					if(ms[j].container!=i) continue;
					innerindex++;
					ibs[innerindex]=new PolygonTypeInnerBoundaryIs(cts[ms[j].index]);
				}
				pgs[polyindex].innerBoundaryIs = ibs;
				ibs=null;
			}

			if(polycount==1)
			{
				fgeo.Add(new gml.PolygonPropertyType(pgs[0]));
			}
			else 
			{
				fgeo.Add(new gml.MultiPolygonPropertyType(new MultiPolygonType(defaultsrs,pgs)));
			}

			ms  = null;
			pgs = null;
			_partPoints = null;
			boxes = null;
		}

		private double[] ConvertToTriangles(double[] orig,int TriType)
		{
			int newcount=3*(orig.Length-2);
			if (newcount<6) return orig;
			if(TriType==(int)PATCHTYPE.TRISTRIP)
			{	
				double[] ts = new double[newcount];
				for(int i=0;i<2;i++) ts[i]=orig[i];	
				for(int i=3;i<orig.Length;i++)
				{
					ts[3*(i-2)]=ts[3*(i-3)+1];	
					ts[3*(i-2)+1]=ts[3*(i-3)+2];
					ts[3*(i-2)+2]=orig[i];
				}
				return ts;
			}
			else if(TriType==(int)PATCHTYPE.TRIFAN)
			{
				double[] ts = new double[newcount];
				for(int i=0;i<2;i++) ts[i]=orig[i];	
				for(int i=3;i<orig.Length;i++)
				{
					ts[3*(i-2)]=ts[3*(i-3)];
					ts[3*(i-2)+1]=ts[3*(i-3)+2];
					ts[3*(i-2)+2]=orig[i];
				}
				return ts;
			}
			return orig;
		}
		#endregion

		#region utilities
		private int SwapEndian32(int orig) 
		{
			//we may also use Encoding on the reader
			byte[] temp = BitConverter.GetBytes(orig);
			Array.Reverse(temp); 
			return BitConverter.ToInt32(temp, 0);
		}

		/// <summary>
		/// Read consecutive integer values
		/// </summary>
		/// <param name="rd">binary reader</param>
		/// <param name="count">number of integer to read</param>
		/// <returns>array of integer value</returns>
		private int[] readIntegerArray(System.IO.BinaryReader rd, int count)
		{
			int bytescount = 4*count;
			int[] da = new int[count];
			byte[] ba = rd.ReadBytes(bytescount);
			for (int n=0;n<da.Length;n++)
			{
				da[n]=BitConverter.ToInt32(ba,n*4);
			}
			ba=null;
			return da;
		}

		/// <summary>
		/// Reads consecutive double values
		/// </summary>
		/// <param name="rd">Binary reader</param>
		/// <param name="count">number of values to read</param>
		/// <returns>array of double value</returns>
		private double[] ReadDoubles(System.IO.BinaryReader rd, int count)
		{
			byte[] ba = rd.ReadBytes(8*count);
			double[] da = new double[count];
			for (int n=0;n<da.Length;n++)
			{
				da[n]=BitConverter.ToDouble(ba,n*8);
			}
			ba=null;
			return da;
		}

		/// <summary>
		/// utility for reading consecutive tuple values
		/// </summary>
		/// <param name="rd">wether we are reading shp or shx</param>
		/// <param name="count">number of tuples</param>
		/// <returns>array of tuples [x0,y0,..,xn,yn]</returns>
		private double[] ReadTuples(System.IO.BinaryReader rd, int tuplesnum)
		{
			return ReadDoubles(rd,tuplesnum*2);
		}
		/// <summary>
		/// utility for reading consecutive tuple values and create the bounding box.
		/// </summary>
		/// <param name="rd">wether we are reading shp or shx</param>
		/// <param name="count">number of tuples</param>
		/// <returns>array of tuples [x0,y0,..,xn,yn]</returns>
		private double[] ReadTuplesWithBox(System.IO.BinaryReader rd, int tuplesnum,ref gml.BoxType box)
		{
			box.Left  = double.MaxValue;
			box.Right = -double.MaxValue;
			box.Bottom = double.MaxValue;
			box.Top = -double.MaxValue;

			byte[] ba = rd.ReadBytes(8*tuplesnum*2);
			double[] da = new double[tuplesnum*2];
			bool odd=true;
			for (int n=0;n<da.Length;n++)
			{
				da[n]=BitConverter.ToDouble(ba,n*8);
				if(odd)
				{
					if(da[n]<box.Left) {box.Left=da[n];}
					else if(da[n]>box.Right) {box.Right=da[n];}
				}
				else
				{
					if(da[n]<box.Bottom) {box.Bottom=da[n];}
					else if(da[n]>box.Top) {box.Top=da[n];}
				}
				odd=!odd;
			}
			ba=null;
			return da;
		}

		private void UpdateStatus(string message, int progressvalue)
		{
			if(progressvalue>0 && progressvalue-StatusArgs.statusValue<5) return;
			StatusArgs.setCurrentStatus(message, progressvalue);
			StatusChange(this,StatusArgs);
		}
		#endregion

	}
}