AE开发第7天

首先先改一个之前的bug，之前做地图代数的时候，用幂函数的式子有错误，应该用Pow而不是Power，很奇怪在arcmap里是Power但是在字符串里只能写Pow，否则会报错。

 private IRaster MapAlgebraOperationSILT(IRasterDataset rasterDataset1, IRasterDataset rasterDataset2)
        {
            IMapAlgebraOp pMapAlgebraOp = new RasterMapAlgebraOpClass();
            IRasterAnalysisEnvironment pEnv = (IRasterAnalysisEnvironment)pMapAlgebraOp;
            pEnv.OutWorkspace = new RasterWorkspaceFactory().OpenFromFile(System.IO.Path.GetDirectoryName(Environment.GetFolderPath(Environment.SpecialFolder.MyPictures)), 0);

            pMapAlgebraOp.BindRaster((IGeoDataset)rasterDataset1, "Ras01");
            pMapAlgebraOp.BindRaster((IGeoDataset)rasterDataset2, "Ras02");
            string sOut = "Pow(([Ras01] / ([Ras01] + [Ras02])), 0.3)";
            return (IRaster)pMapAlgebraOp.Execute(sOut);
        }

就是这一段代码报错。

今天做的是降水侵蚀度的计算，用到了地图代数和重分类的内容。

 private IRaster reclass2(IRaster pRaster)
        {
            IRasterProps rasterProps = (IRasterProps)pRaster;

            //设置栅格数据起始点  
            IPnt pBlockSize = new Pnt();
            pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);

            //选取整个范围  
            IPixelBlock pPixelBlock = pRaster.CreatePixelBlock(pBlockSize);

            //左上点坐标  
            IPnt tlp = new Pnt();
            tlp.SetCoords(0, 0);

            //读入栅格  
            IRasterBandCollection pRasterBands = pRaster as IRasterBandCollection;
            IRasterBand pRasterBand = pRasterBands.Item(0);
            IRawPixels pRawPixels = pRasterBands.Item(0) as IRawPixels;
            pRawPixels.Read(tlp, pPixelBlock);

            //将PixBlock的值组成数组  
            System.Array pSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;
            List<double> pixelValues = new List<double>();
            System.Array dataArray = pSafeArray;

            for (int y = 0; y < rasterProps.Height; y++)
            {
                for (int x = 0; x < rasterProps.Width; x++)
                {
                    double value = Convert.ToDouble(dataArray.GetValue(x, y));

                    if (value <= 0)
                        dataArray.SetValue(0.0f, x, y);
                    else if (value > 0)
                    {
                        dataArray.SetValue((float)value, x, y);
                    }
                    else
                    {
                        dataArray.SetValue(float.NaN, x, y); // 设置为空值
                    }
                }
            }
            pPixelBlock.set_SafeArray(0, pSafeArray);

            //编辑raster，将更新的值写入raster中  
            IRasterEdit rasterEdit = pRaster as IRasterEdit;
            rasterEdit.Write(tlp, pPixelBlock);
            rasterEdit.Refresh();
            return pRaster;
        }
        private void r因子计算ToolStripMenuItem1_Click(object sender, EventArgs e)
        {
            // 自定义的提示信息
            string customMessage1 = "请输入年平均降水栅格";
            // 显示自定义提示信息
            MessageBox.Show(customMessage1, "提示", MessageBoxButtons.OK, MessageBoxIcon.Information);
            IRasterDataset rasterDataset = OpenTiffFile();
            if (rasterDataset == null)
            {
                MessageBox.Show("未选择任何TIFF文件或者打开文件失败。");
                return;
            }
            // 自定义的提示信息
            string customMessage2 = "请输入月平均降水栅格";
            // 显示自定义提示信息
            MessageBox.Show(customMessage2, "提示", MessageBoxButtons.OK, MessageBoxIcon.Information);
            // 用户选择的多个TIFF文件
            List<IRasterDataset> rasterDatasets = OpenTiffFiles();
            if (rasterDatasets == null)
            {
                MessageBox.Show("未选择任何TIFF文件或者打开文件失败。");
                return;
            }
            IRaster resultRaster = RR(rasterDatasets, rasterDataset);
            IRasterLayer rasterLayer = new RasterLayerClass();
            rasterLayer.CreateFromRaster(resultRaster);
            rasterLayer.Name = "R因子计算";

            axMapControl2.AddLayer((ILayer)rasterLayer, 0);
            axMapControl2.ActiveView.Refresh();
            SaveRasterToFile(resultRaster);
        }
        private IRaster RR(List<IRasterDataset> rasterDatasets, IRasterDataset rasterDataset)
        {
            if (rasterDatasets.Count == 0 || rasterDataset == null)
            {
                return null;
            }

            // 获取第一个栅格数据集作为基础进行相加
            IRasterDataset baseDataset = rasterDataset;
            IMapAlgebraOp pMapAlgebraOp = new RasterMapAlgebraOpClass();
            IRasterAnalysisEnvironment pEnv = (IRasterAnalysisEnvironment)pMapAlgebraOp;
            pEnv.OutWorkspace = new RasterWorkspaceFactory().OpenFromFile(System.IO.Path.GetDirectoryName(Environment.GetFolderPath(Environment.SpecialFolder.MyPictures)), 0);
            IGeoDataset baseDataset1 = baseDataset as IGeoDataset;

            // 创建结果栅格数据集列表
            List<IRasterDataset> rasterDatasets1 = new List<IRasterDataset>();
            // 绑定栅格数据
            pMapAlgebraOp.BindRaster(baseDataset1, "Ras01");
            for (int i = 0; i < rasterDatasets.Count; i++)
            {
                IRasterDataset currentDataset = rasterDatasets[i];
                pMapAlgebraOp.BindRaster((IGeoDataset)currentDataset, "Ras02");

                // 设置地图代数表达式
                string sOut = "1.5 * Log10(([Ras02] + 0.001) * ([Ras02] + 0.001) / [Ras01])";

                // 执行栅格操作
                IGeoDataset result = pMapAlgebraOp.Execute(sOut);
                IRaster rasterResult = (IRaster)result;
                IRaster rasterResult1 = reclass2(rasterResult);
                // 将 IRaster 转换为 IRasterDataset
                IRasterDataset rasterDatasetResult = (rasterResult1 as IRaster2).RasterDataset;

                // 将结果更新为新的基础栅格数据集，以便进行下一轮相加
                rasterDatasets1.Add(rasterDatasetResult);
            }

            IGeoDataset newDataset1 = rasterDatasets1[0] as IGeoDataset;

            // 从第二个栅格数据集开始进行相加
            for (int i = 1; i < rasterDatasets1.Count; i++)
            {
                // 绑定栅格数据
                pMapAlgebraOp.BindRaster(newDataset1, "Ras01");
                IRasterDataset currentDataset = rasterDatasets1[i];
                pMapAlgebraOp.BindRaster((IGeoDataset)currentDataset, "Ras02");

                // 设置地图代数表达式
                string sOut = "[Ras01] + [Ras02]";

                // 创建栅格运算对象
                IGeoDataset result = pMapAlgebraOp.Execute(sOut);

                // 如果栅格相加操作失败，则返回null
                if (result == null)
                {
                    return null;
                }

                // 将结果更新为新的基础栅格数据集，以便进行下一轮相加
                newDataset1 = result;
            }
            return (IRaster)newDataset1; // 返回相加后的栅格数据集
        }

最后做出来的效果如下