/// <summary> /// 线与TIN相交,获取剖面样点的x,y,z /// </summary> /// <param name="pPolyLine">剖面线</param> /// <param name="pSurface">Tin表面</param> /// <param name="interpolatePtsCount">生成多少个插值点</param> public IPointCollection RunTinProfile(IPolyline pPolyLine, ITinSurface pSurface,int interpolatePtsCount) { IEnvelope pSurfaceEnvelope, pPolyLineEnvelope;//Envelop相当于要素的范围,东西南北 pPolyLineEnvelope = new EnvelopeClass(); pPolyLine.QueryEnvelope(pPolyLineEnvelope); pSurfaceEnvelope = pSurface.Domain.Envelope; IPointCollection pMultiPoint = new MultipointClass(); if (pSurfaceEnvelope.XMax < pPolyLineEnvelope.XMax || pSurfaceEnvelope.XMin > pPolyLineEnvelope.XMin || pSurfaceEnvelope.YMax < pPolyLineEnvelope.YMax || pSurfaceEnvelope.YMin > pPolyLineEnvelope.YMin) { MessageBox.Show("The Polyline you drew (selected) is out of surface layer's(the active layer) boundary!!! "); return pMultiPoint; } IPoint outPoint; IZAware pZA, pZAM; int i; double thePstn; bool bAsRatio; ICurve pCurve = pPolyLine as ICurve; ITinSurfaceElement pTinSurfaceElement; pZAM = pMultiPoint as IZAware; pZAM.ZAware = true; object missing = Type.Missing; for (i = 0; i <= interpolatePtsCount - 1; i++) { outPoint = new PointClass(); thePstn = ((double)i / (double)(interpolatePtsCount - 1)); //参数2和参数3配套使用, 参数2代表距离,参数3指示2是否按Curve的长度比率来 pCurve.QueryPoint(esriSegmentExtension.esriNoExtension, thePstn, true, outPoint); pTinSurfaceElement = pSurface.GetSurfaceElement(outPoint); pZA = outPoint as IZAware; pZA.ZAware = true; outPoint.Z = pTinSurfaceElement.Elevation; pMultiPoint.AddPoint(outPoint, ref missing, ref missing);//点加在哪个点之前,哪个点之后 } return pMultiPoint; }
/// <summary> /// 剖面线 /// </summary> /// <param name="sender"></param> /// <param name="e"></param> private void butn_profile_Click(object sender, EventArgs e) { ITin pTin = GetTinLayer(@"E:\等高线\", "tin"); Tin tempTin=new TinClass(); RunTinProfile(this._profileline, pTin as ITinSurface, 100); }
//基于Raster
/// <summary> /// 剖面分析 /// 基于Raster数据 ///PointArray的点是真实坐标,不是屏幕坐标 ///备注：尚未测试多个Raster的投影面积不同的情况 /// </summary> /// <param name="pPolyLine">剖面线</param> /// <param name="pRaster">raster数据</param> /// <param name="interpolatePtsCount">插值点数</param> /// <param name="lineLength">剖面线长度 如果只是部分相交,那么会被更新</param> /// <param name="partialintersect">是否只部分相交</param> /// <param name="pPartPolyline">相交的剖线</param> /// <returns></returns> public static IPointArray RunRasterProfile(IPolyline pPolyLine, IRaster pRaster, int interpolatePtsCount, ref double lineLength, out bool partialintersect, out IPolyline pPartPolyline) { partialintersect = false; pPartPolyline = new PolylineClass();//相交部分的线 IPointArray pPointArray = new PointArrayClass(); IEnvelope pRasterEnvelope, pPolyLineEnvelope;//Envelop相当于要素的范围,东西南北 pPolyLineEnvelope = new EnvelopeClass(); pRasterEnvelope=GetEnvelopeOfRaster(pRaster); IRaster2 pRaster2 = pRaster as IRaster2; IPointCollection pMultiPoint = new MultipointClass(); if (pRasterEnvelope.XMax < pPolyLine.Envelope.XMax || pRasterEnvelope.XMin > pPolyLine.Envelope.XMin || pRasterEnvelope.YMax < pPolyLine.Envelope.YMax || pRasterEnvelope.YMin > pPolyLine.Envelope.YMin) { pPolyLine = GetClippedPolyline(pPolyLine, pRasterEnvelope);//获取相交部分的Polyline if (pPolyLine.Length == 0)//剖线没有经过任何钻孔数据区 { return pPointArray; } else //只是部分相交 { pPartPolyline = pPolyLine;//相交部分的线 lineLength = pPolyLine.Length; partialintersect = true; } } else pPartPolyline = pPolyLine;//完全在范围内 IPoint outPoint; int i; double thePstn; ICurve pCurve = pPolyLine as ICurve; object missing = Type.Missing; for (i = 0; i <= interpolatePtsCount - 1; i++) { try { outPoint = new PointClass(); thePstn = ((double)i / (double)(interpolatePtsCount - 1)); //参数2和参数3配套使用, 参数2代表距离,参数3指示2是否按Curve的长度比率来 //获取表面上某点的坐标 //pCurve.QueryPoint(esriSegmentExtension.esriExtendTangentAtFrom, thePstn, true, outPoint); pCurve.QueryPoint(esriSegmentExtension.esriNoExtension, thePstn, true, outPoint); double pixelvalue = IdentifyPixelValue(pRaster, outPoint.X, outPoint.Y);//根据坐标获取 该点的pixelvalue if (pixelvalue == -99999) continue; //如果改点pixel没有值 outPoint.Z = pixelvalue; ESRI.ArcGIS.Geometry.Point pPoint = new ESRI.ArcGIS.Geometry.Point(); pPoint.X = thePstn * lineLength; pPoint.Z = outPoint.Z;//点的高程 pPointArray.Add(pPoint as IPoint); } catch { continue; } } return pPointArray; }
/// <summary> /// 从文件获取Raster /// </summary> /// <param name="folder">Raster的父目录</param> /// <param name="tinDataset">Raster目录</param> /// <returns></returns> public static IRaster GetRaster(string folder, string rasterName) { IWorkspaceFactory workspaceFactory = new RasterWorkspaceFactory(); IWorkspace workspace; workspace = workspaceFactory.OpenFromFile(folder, 0); //inPath栅格数据存储路径 IRasterWorkspace pRasterWS = workspace as IRasterWorkspace; IRasterDataset pRasterDataSet= pRasterWS.OpenRasterDataset(rasterName); IRaster pRaster= pRasterDataSet.CreateDefaultRaster(); return pRaster; }
/// <summary> /// 获取Raster的Envelope /// </summary> /// <returns></returns> public static IEnvelope GetEnvelopeOfRaster(IRaster pRaster) { IRasterProps rasterProps = (IRasterProps)pRaster; IEnvelope extent = rasterProps.Extent; return extent; }
/// <summary> /// 根据坐标获取栅格数据的PixelValue /// </summary> /// <param name="raster"></param> /// <param name="xMap"></param> /// <param name="yMap"></param> public static double IdentifyPixelValue(IRaster raster, double xMap, double yMap) { IRaster2 raster2 = (IRaster2)raster; int row = 0; int col = 0; raster2.MapToPixel(xMap, yMap,out col,out row); double pixelValue; //Get the value at a given band. if (raster2.GetPixelValue(0, col, row) == null) pixelValue = -99999; else pixelValue = Convert.ToDouble(raster2.GetPixelValue(0, col, row)); /*Convert.ToDouble(null)=0.0 * 但GetPixelValue的返回确实可能为null */ return pixelValue; }
/// <summary> /// rectangle/envelop 裁剪线 /// </summary> /// <param name="pPolyLine">剖面线</param> /// <param name="pEnvelop">Tin的Envelop</param> /// <returns></returns> public static IPolyline GetClippedPolyline(IPolyline pPolyLine,IEnvelope pEnvelop) { ISpatialReferenceFactory spatialReferenceFactory = new SpatialReferenceEnvironmentClass(); //如果polyline和Envelop是GCS的 //IGeographicCoordinateSystem gcsSys = spatialReferenceFactory.CreateGeographicCoordinateSystem((int)esriSRGeoCS3Type.esriSRGeoCS_Xian1980); //pPolyLine.SpatialReference = gcsSys; //pEnvelop.SpatialReference = gcsSys; //如果polyline和Envelop是PCS的 IProjectedCoordinateSystem pcsSys = spatialReferenceFactory.CreateProjectedCoordinateSystem((int)esriSRProjCS4Type.esriSRProjCS_Xian1980_3_Degree_GK_Zone_39); pPolyLine.SpatialReference = pcsSys; pEnvelop.SpatialReference = pcsSys; ITopologicalOperator2 pTopOper = pPolyLine as ITopologicalOperator2; IGeometry pGeoafter=new PolylineClass(); pTopOper.QueryClipped(pEnvelop, pGeoafter); return pGeoafter as IPolyline; }