首先先改一个之前的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; // 返回相加后的栅格数据集
}
最后做出来的效果如下