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

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 "";
}
}