进行波段运算_C#结合GDAL实现对遥感影像进行波段运算

一、编译安装GDAL

按如下教程编译GDAL，所选的IDE为vs2010，编译时正确：

https://blog.csdn.net/zhangying_496/article/details/51258862

在使用开发C#下的GDAL时，编译时正确。但在运行程序是出现了如下问题：

Gdal.AllRegister();

出现异常：OSGeo.GDAL.GdalPINVOKE”的类型初始值设定项引发异常。

解决方法：编译正常，说明已经正确添加GDAL主要的八个动态库DLL：gdal_csharp.dll，gdal_wrap.dll，gdalconst_csharp，gdalconst_wrap.dll，ogr_csharp.dll，ogr_wrap.dll，osr_csharp.dll，osr_wrap.dll。

最后从以下网站下载了编译好的GDAL

GDAL C# DLL下载：

http://www.gisinternals.com/sdk/​www.gisinternals.com

在这里我下载了release-1600-x64-gdal-1-10-1-mapserver-6-4-1

解压后，将bin文件夹里的N个dll，文件拷贝在执行debug目录下即可解决上面遇到的问题。

二、读取遥感影像、波段计算代码

OSGeo.GDAL.Dataset ds;
 public void button1_Click(object sender, EventArgs e)
        {
            OSGeo.GDAL.Gdal.AllRegister();
            ds = OSGeo.GDAL.Gdal.Open(textBox1.Text, OSGeo.GDAL.Access.GA_Update);
 int srcWidth = ds.RasterXSize;
 int srcHeight = ds.RasterYSize;
 MessageBox.Show("读取影像成功");
 int m = Convert.ToInt16(textBox2.Text);
 int n = Convert.ToInt16(textBox3.Text);
            OSGeo.GDAL.Band srcBand1 = ds.GetRasterBand(m);
            OSGeo.GDAL.Band srcBand2 = ds.GetRasterBand(n);
            OSGeo.GDAL.DataType srcType = ds.GetRasterBand(m).DataType;
 double val;
 int hasval;
            srcBand1.GetMaximum(out val,out hasval);
 double[] tmpD = new double[6];
            ds.GetGeoTransform(tmpD);
 double[] newArray = new double[srcWidth * srcHeight];
 double[] dataArray1 = new double[srcWidth * srcHeight];
 double[] dataArray2 = new double[srcWidth * srcHeight];
 int nXSize = ds.RasterXSize;
            srcBand1.ReadRaster (0, 0, srcWidth, srcHeight, dataArray1, srcWidth, srcHeight, 0, 0);
            srcBand2.ReadRaster(0, 0, srcWidth, srcHeight, dataArray2, srcWidth, srcHeight, 0, 0);
            OSGeo.GDAL.Driver dstDrv = Gdal.GetDriverByName("GTiff");
 int counts=srcWidth * srcHeight;
 for(int icount = 0 ; icount<counts ; icount++)
            {
                newArray[icount] = (dataArray1[icount] - dataArray2[icount]) / (dataArray1[icount] + dataArray2[icount]+1e-7);
             }
 double yz = Convert.ToDouble(textBox5.Text);
 for (int icount = 0; icount < counts; icount++)
            {
 if (newArray[icount] > yz)
                    newArray[icount] = 1;
 else
                    newArray[icount] = 0;
            }
            OSGeo.GDAL.Dataset multiDstDataset = dstDrv.Create(textBox4.Text, srcWidth, srcHeight, 1, OSGeo.GDAL.DataType.GDT_Float32, null);
            multiDstDataset.SetGeoTransform(tmpD);
            multiDstDataset.SetProjection(ds.GetProjection()); 
 int multiBandCount = 1;
 int[] bandNum = new int[multiBandCount];
 for (int i = 0; i < multiBandCount; ++i)
            {
                bandNum[i] = i + 1;
            }
                       multiDstDataset.WriteRaster(0, 0, srcWidth, srcHeight, newArray, srcWidth, srcHeight, 1, bandNum, 0, 0, 0);
            multiDstDataset.FlushCache();
 MessageBox.Show("影像已成功保存");
        }

处理得到结果

下一步计划是将二值影像转换为shapefile格式。