AE 实现3DAnalysis Profile Graph功能 （基于TIN,基于Raster 得到线上所有点的高程）...

  /// <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;
        }