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

来源：互联网发布：网易mumu模拟器mac版编辑：程序博客网时间：2024/06/05 18:10
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 "";            }        }
0 0