AE + GDAL实现影像按标准图幅分割(下)

时间:2023-03-08 15:40:36

  在上篇实现了遥感影像的切割,本篇讲切割前的准备。主要分为以下几步:

  (1)将影像的投影坐标转换为地理坐标,以便于之后的图幅划分。AE坐标转换函数如下

 private bool Proj2Geo(ISpatialReference pspr, double xProj, double yProj, ref double xGeo, ref double yGeo)
{
if (pspr is IGeographicCoordinateSystem)
return false;
IProjectedCoordinateSystem pcs = pspr as IProjectedCoordinateSystem;
ISpatialReference gspr = pcs.GeographicCoordinateSystem; IPoint pt = new PointClass();
pt.PutCoords(xProj, yProj);
pt.SpatialReference = pspr;
pt.Project(gspr);
xGeo = pt.X;
yGeo = pt.Y; return true;
} private bool Geo2Proj(ISpatialReference pspr, double xGeo, double yGeo, ref double xProj, ref double yProj)
{
if (pspr is IGeographicCoordinateSystem)
return false;
IProjectedCoordinateSystem pcs = pspr as IProjectedCoordinateSystem;
ISpatialReference gspr = pcs.GeographicCoordinateSystem; IPoint pt = new PointClass();
pt.PutCoords(xGeo, yGeo);
pt.SpatialReference = gspr;
pt.Project(pspr);
xProj = pt.X;
yProj = pt.Y; return true;
}

  (2)计算包含影像的标准图幅的四至,首先我定义了一个格网类

 public class GridInfo
{
public double XMin;
public double XMax;
public double YMin;
public double YMax;
public int rows;
public int cols; public void setGridInfo(double xMax, double xMin, double yMax, double yMin, int rowCount, int colCount)
{
XMax = xMax;
XMin = xMin;
YMax = yMax;
YMin = yMin;
rows = rowCount;
cols = colCount;
}
}

根据比例尺计算格网的大小

 private GridInfo setGridInfoByScale(int scale)
{
GridInfo gridInfo = new GridInfo();
double dxInSnd, dyInSnd;
switch (scale)
{
case :
dxInSnd = Degree.Degree2Second(0.03125);
dyInSnd = Degree.Minute2Second(1.25);
break;
case :
dxInSnd = Degree.Degree2Second(0.0625);
dyInSnd = Degree.Minute2Second(2.5);
break;
case :
dxInSnd = Degree.Minute2Second(7.5);
dyInSnd = Degree.Minute2Second();
break;
case :
dxInSnd = Degree.Minute2Second();
dyInSnd = Degree.Minute2Second();
break;
case :
dxInSnd = Degree.Minute2Second();
dyInSnd = Degree.Minute2Second();
break;
case :
dxInSnd = Degree.Degree2Second(1.5);
dyInSnd = Degree.Degree2Second();
break;
case :
dxInSnd = Degree.Degree2Second();
dyInSnd = Degree.Degree2Second();
break;
case :
dxInSnd = Degree.Degree2Second();
dyInSnd = Degree.Degree2Second();
break;
default:
dxInSnd = ;
dyInSnd = ;
break;
} if (dxInSnd == dyInSnd && dxInSnd == 0.0)
return null; double pXMin = 0.0, pXMax = 0.0, pYMin = 0.0, pYMax = 0.0;
double gXMin = 0.0, gXMax = 0.0, gYMin = 0.0, gYMax = 0.0;
if (envelope.SpatialReference is IProjectedCoordinateSystem)
{
Proj2Geo(envelope.SpatialReference, envelope.XMax, envelope.YMax, ref gXMax, ref gYMax);
Proj2Geo(envelope.SpatialReference, envelope.XMin, envelope.YMin, ref gXMin, ref gYMin);
}
else
{
gXMax = envelope.XMax; gXMin = envelope.XMin;
gYMax = envelope.YMax; gYMin = envelope.YMin;
}
int cols =Convert.ToInt32(Math.Round(Degree.Degree2Second( - gXMin) / dxInSnd)) + ;
gridInfo.XMin = - Degree.Second2Degree(cols * dxInSnd);
cols = Convert.ToInt32(Math.Round(Degree.Degree2Second( - gXMax) / dxInSnd));
gridInfo.XMax = - Degree.Second2Degree(cols * dxInSnd); int rows = Convert.ToInt32(Math.Round(Degree.Degree2Second(gYMin) / dyInSnd));
gridInfo.YMin =Degree.Second2Degree(rows * dyInSnd);
rows = Convert.ToInt32(Math.Round(Degree.Degree2Second(gYMax) / dyInSnd)) + ;
gridInfo.YMax =Degree.Second2Degree( rows * dyInSnd); gridInfo.rows = Convert.ToInt32(Math.Round(Degree.Degree2Second(gridInfo.YMax - gridInfo.YMin) / dyInSnd));
gridInfo.cols = Convert.ToInt32(Math.Round(Degree.Degree2Second(gridInfo.XMax - gridInfo.XMin) / dxInSnd)); if (envelope.SpatialReference is IProjectedCoordinateSystem)
{
Geo2Proj(envelope.SpatialReference, gridInfo.XMax, gridInfo.YMax, ref pXMax, ref pYMax);
Geo2Proj(envelope.SpatialReference, gridInfo.XMin, gridInfo.YMin, ref pXMin, ref pYMin);
gridInfo.XMax = pXMax; gridInfo.XMin = pXMin;
gridInfo.YMax = pYMax; gridInfo.YMin = pYMin;
} return gridInfo;
}

  (3)图幅命名类

 public class FrameName
{
private int scale;
/// <summary>
///
/// </summary>
/// <param name="Scale">分数比例尺的分母(比例尺为1:50000则参数为50000)</param>
public FrameName(int Scale)
{
scale = Scale;
} public String getFrameName(double longtitude, double latitude)
{
int baseRowCount=(int)(latitude / );
char baseChar = Convert.ToChar('A' + baseRowCount);
int baseNum = (int)(longtitude / ) + ; if (scale == )
return String.Format("{0}{1}", baseChar, baseNum); char scaleChar = 'E'; //1:50000比例尺的代码为E,其他以后补充
switch (scale)
{
case :
scaleChar='B';break;
case :
scaleChar = 'C';break;
case :
scaleChar = 'D';break;
case :
scaleChar='E';break;
case :
scaleChar = 'F';break;
case :
scaleChar = 'G';break;
case :
scaleChar = 'H';break;
} double dxInSnd, dyInSnd;
switch (scale)
{
case :
dxInSnd = 0.03125;
dyInSnd = 0.0208333333;
break;
case :
dxInSnd = 0.0625;
dyInSnd = 0.0416666667;
break;
case :
dxInSnd = 0.125;
dyInSnd = 0.0833333333;
break;
case :
dxInSnd = 0.25;
dyInSnd = 0.1666666667;
break;
case :
dxInSnd = 0.5;
dyInSnd = 0.3333333333;
break;
case :
dxInSnd = 1.5;
dyInSnd = ;
break;
case :
dxInSnd = ;
dyInSnd = ;
break;
default:
dxInSnd = ;
dyInSnd = ;
break;
} if (dxInSnd == dyInSnd && dxInSnd == 0.0)
return null; int row = (int)(((baseRowCount + ) * - latitude)/dyInSnd) + ;
//int row = 24-((int)(latitude / dyInSnd)) % 24;
int col = (int)((longtitude % ) /dxInSnd) + ; return String.Format("{0}{1}{2}{3}{4}", baseChar, baseNum, scaleChar, row.ToString().PadLeft(, ''), col.ToString().PadLeft(, ''));
} }

  (4)调用函数进行分割

 private void CreateFrame(GridInfo gridInfo, String outDir, int scale)
{
double xLength = (gridInfo.XMax - gridInfo.XMin) / gridInfo.cols;
double yLength = (gridInfo.YMax - gridInfo.YMin) / gridInfo.rows;
FrameName frameName = new FrameName(scale); for (int i = ; i < gridInfo.rows; i++)
{
double yMin = gridInfo.YMin + i * yLength;
double yMax = gridInfo.YMin + (i + ) * yLength;
for (int j = ; j < gridInfo.cols; j++)
{
double xMin = gridInfo.XMin + j * xLength;
double xMax = gridInfo.XMin + (j + ) * xLength;
String pszOutFile; if (envelope.SpatialReference is IProjectedCoordinateSystem)
{
double prjX = 0.0, prjY = 0.0;
Proj2Geo(envelope.SpatialReference, (xMax + xMin) / , (yMax + yMin) / , ref prjX, ref prjY);
pszOutFile = frameName.getFrameName(prjX,prjY);
}
else
pszOutFile = frameName.getFrameName((xMax + xMin) / , (yMax + yMin) / );
pszOutFile = outDir + "\\" + pszOutFile + ".tif";
CutUpdateImageByAOIGDAL(rasterLyr.FilePath, pszOutFile, xMin, xMax, yMin, yMax, "GTiff");
}
}
}

  通过这4步,一个简易的图幅分割工具就做好了。