C#应用GDAL通过传入范围获取tif数据的最大高程值及其坐标

GDAL的安装编译可参考 之前的博文。点击打开链接

此函数可通过传入一个规则范围position=”left,top,rigth,bottom”,返回这个范围内的最大高程及其坐标和最小高程及其坐标

 public string GetMultifyElevation(string positions)
        {
            positions = "116.0,40.166667,116.25,40.0";//模拟传入的范围（left，top，right，bottom）
            float dProjX, dProjY , dProjX1, dProjY1;
            string[] arr = positions.Split(',');
            dProjX =float.Parse(arr[0]);
            dProjY = float.Parse(arr[1]);
            dProjX1 = float.Parse(arr[2]);
            dProjY1 = float.Parse(arr[3]);


            string strFilePath = "C:\\webservices\\data\\srtm\\chinaclip.tif";//读取的文件路径
            string testPath = "C:\\webservices\\data\\chinaclip.tif";//要写的文件路径
            string elvate = "";
            Gdal.AllRegister();
            Gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");    //支持中文
            Dataset ds = Gdal.Open(strFilePath, Access.GA_ReadOnly);
            try
            {
                Band Band = ds.GetRasterBand(1);


                //获取图像的尺寸               
                int width = Band.XSize;
                int height = Band.YSize;




                //获取坐标变换系数
                double[] adfGeoTransform = new double[6];
                ds.GetGeoTransform(adfGeoTransform);
                //获取行列号
                double dTemp = adfGeoTransform[1] * adfGeoTransform[5] - adfGeoTransform[2] * adfGeoTransform[4];
                double dCol = 0.0, dRow = 0.0;
              


                    dCol = (adfGeoTransform[5] * (dProjX - adfGeoTransform[0]) -
                        adfGeoTransform[2] * (dProjY - adfGeoTransform[3])) / dTemp + 0.5;
                    dRow = (adfGeoTransform[1] * (dProjY - adfGeoTransform[3]) -
                        adfGeoTransform[4] * (dProjX - adfGeoTransform[0])) / dTemp + 0.5;
                    int dc = Convert.ToInt32(dCol);
                    int dr = Convert.ToInt32(dRow);


                    dCol = (adfGeoTransform[5] * (dProjX1 - adfGeoTransform[0]) -
                           adfGeoTransform[2] * (dProjY1 - adfGeoTransform[3])) / dTemp + 0.5;
                    dRow = (adfGeoTransform[1] * (dProjY1 - adfGeoTransform[3]) -
                        adfGeoTransform[4] * (dProjX1 - adfGeoTransform[0])) / dTemp + 0.5;
                    int dc1 = Convert.ToInt32(dCol);
                    int dr1 = Convert.ToInt32(dRow);
                    int fx = dc - dc1;
                    int fy = dr - dr1;
                    if (fx < 0)
                        fx = -fx;
                    if (fy < 0)
                        fy = -fy;


                    //获取DEM数值到一维数组
                    float[] data = new float[fx * fy];
                    //读取感觉可以去掉
                    CPLErr err = Band.ReadRaster(dc, dr, fx, fy, data, fx, fy, 0, 0);
                    //裁切
                    int[] bandmap = { 1 };
                    DataType DT = DataType.GDT_CFloat32;
                    Dataset dataset = ds.GetDriver().Create(testPath, fx, fy, 1, DT, null);
                    //写入
                    dataset.WriteRaster(0, 0, fx, fy, data, fx, fy, 1, bandmap, 0, 0, 0);


                   


                    Band bd = dataset.GetRasterBand(1);
                    //获取最小最大值
                    double[] lim = new double[2];
                    bd.ComputeRasterMinMax(lim, 1);   //最后的ApproxOK设为1，设为0效果一样
                    float _Min = (float)lim[0];
                    float _Max = (float)lim[1];   


                    bd.ReadRaster(0, 0, fx, fy, data, fx, fy, 0, 0); 
                int c =0;
                int x1 = -1, y1 = -1, x2 = -1, y2 = -1;
                //遍历获取行列值
                    for (int i = 0; i < fx; i++)
                    {
                        if (x2 != -1 && y2 != -1 && x1 != -1 && y1 != -1)
                        { 
                            goto conmehere;//为了提高效率使用goto语句
                        }


                        for (int j = 0; j < fy; j++)
                        {
                            if (x2 != -1 && y2 != -1 && x1 != -1 && y1 != -1)
                            {
                                goto conmehere;
                            }
                            if (_Min == data[c++])
                            {


                                x1 = i;
                                y1 = j;                               
                            }
                            if (_Max == data[c])
                            {
                                x2 = i;
                                y2 = j;                               
                            }
                        }
                    }
                conmehere:


                    //adfGeoTransform[0]  左上角x坐标   
                    //adfGeoTransform[1]  东西方向分辨率  
                    //adfGeoTransform[2]  旋转角度, 0表示图像 "北方朝上"  
                    //adfGeoTransform[3]  左上角y坐标   
                    //adfGeoTransform[4]  旋转角度, 0表示图像 "北方朝上"  
                    //adfGeoTransform[5]  南北方向分辨率  
                    Band.Dispose();
                    double geox1 = dProjX + x1 * adfGeoTransform[1] + y1 * adfGeoTransform[2];
                    double geoy1 = dProjY + x1 * adfGeoTransform[4] + y1 * adfGeoTransform[5];
                    double geox2 = dProjX + x2 * adfGeoTransform[1] + y2 * adfGeoTransform[2];
                    double geoy2 = dProjY + x2 * adfGeoTransform[4] + y2 * adfGeoTransform[5];


                    elvate = geox1 + "," + geoy1 + "," + _Min + ";" + geox2 + "," + geoy2 + "," + _Max;
                return elvate;
            }
            catch
            {
                return "";
            }
        }

转载自：https://blog.csdn.net/nothing_is_imposible/article/details/18313417