在上篇实现了遥感影像的切割,本篇讲切割前的准备。主要分为以下几步:
(1)将影像的投影坐标转换为地理坐标,以便于之后的图幅划分。AE坐标转换函数如下
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 private bool Proj2Geo(ISpatialReference pspr, double xProj, double yProj, ref double xGeo, ref double yGeo) 2 { 3 if (pspr is IGeographicCoordinateSystem) 4 return false; 5 IProjectedCoordinateSystem pcs = pspr as IProjectedCoordinateSystem; 6 ISpatialReference gspr = pcs.GeographicCoordinateSystem; 7 8 IPoint pt = new PointClass(); 9 pt.PutCoords(xProj, yProj); 10 pt.SpatialReference = pspr; 11 pt.Project(gspr); 12 xGeo = pt.X; 13 yGeo = pt.Y; 14 15 return true; 16 } 17 18 private bool Geo2Proj(ISpatialReference pspr, double xGeo, double yGeo, ref double xProj, ref double yProj) 19 { 20 if (pspr is IGeographicCoordinateSystem) 21 return false; 22 IProjectedCoordinateSystem pcs = pspr as IProjectedCoordinateSystem; 23 ISpatialReference gspr = pcs.GeographicCoordinateSystem; 24 25 IPoint pt = new PointClass(); 26 pt.PutCoords(xGeo, yGeo); 27 pt.SpatialReference = gspr; 28 pt.Project(pspr); 29 xProj = pt.X; 30 yProj = pt.Y; 31 32 return true; 33 }
(2)计算包含影像的标准图幅的四至,首先我定义了一个格网类
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 public class GridInfo 2 { 3 public double XMin; 4 public double XMax; 5 public double YMin; 6 public double YMax; 7 public int rows; 8 public int cols; 9 10 public void setGridInfo(double xMax, double xMin, double yMax, double yMin, int rowCount, int colCount) 11 { 12 XMax = xMax; 13 XMin = xMin; 14 YMax = yMax; 15 YMin = yMin; 16 rows = rowCount; 17 cols = colCount; 18 } 19 }
根据比例尺计算格网的大小
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 private GridInfo setGridInfoByScale(int scale) 2 { 3 GridInfo gridInfo = new GridInfo(); 4 double dxInSnd, dyInSnd; 5 switch (scale) 6 { 7 case 5000: 8 dxInSnd = Degree.Degree2Second(0.03125); 9 dyInSnd = Degree.Minute2Second(1.25); 10 break; 11 case 10000: 12 dxInSnd = Degree.Degree2Second(0.0625); 13 dyInSnd = Degree.Minute2Second(2.5); 14 break; 15 case 25000: 16 dxInSnd = Degree.Minute2Second(7.5); 17 dyInSnd = Degree.Minute2Second(5); 18 break; 19 case 50000: 20 dxInSnd = Degree.Minute2Second(15); 21 dyInSnd = Degree.Minute2Second(10); 22 break; 23 case 100000: 24 dxInSnd = Degree.Minute2Second(30); 25 dyInSnd = Degree.Minute2Second(20); 26 break; 27 case 250000: 28 dxInSnd = Degree.Degree2Second(1.5); 29 dyInSnd = Degree.Degree2Second(1); 30 break; 31 case 500000: 32 dxInSnd = Degree.Degree2Second(3); 33 dyInSnd = Degree.Degree2Second(2); 34 break; 35 case 1000000: 36 dxInSnd = Degree.Degree2Second(6); 37 dyInSnd = Degree.Degree2Second(4); 38 break; 39 default: 40 dxInSnd = 0; 41 dyInSnd = 0; 42 break; 43 } 44 45 if (dxInSnd == dyInSnd && dxInSnd == 0.0) 46 return null; 47 48 double pXMin = 0.0, pXMax = 0.0, pYMin = 0.0, pYMax = 0.0; 49 double gXMin = 0.0, gXMax = 0.0, gYMin = 0.0, gYMax = 0.0; 50 if (envelope.SpatialReference is IProjectedCoordinateSystem) 51 { 52 Proj2Geo(envelope.SpatialReference, envelope.XMax, envelope.YMax, ref gXMax, ref gYMax); 53 Proj2Geo(envelope.SpatialReference, envelope.XMin, envelope.YMin, ref gXMin, ref gYMin); 54 } 55 else 56 { 57 gXMax = envelope.XMax; gXMin = envelope.XMin; 58 gYMax = envelope.YMax; gYMin = envelope.YMin; 59 } 60 int cols =Convert.ToInt32(Math.Round(Degree.Degree2Second(180 - gXMin) / dxInSnd)) + 1; 61 gridInfo.XMin = 180 - Degree.Second2Degree(cols * dxInSnd); 62 cols = Convert.ToInt32(Math.Round(Degree.Degree2Second(180 - gXMax) / dxInSnd)); 63 gridInfo.XMax = 180 - Degree.Second2Degree(cols * dxInSnd); 64 65 int rows = Convert.ToInt32(Math.Round(Degree.Degree2Second(gYMin) / dyInSnd)); 66 gridInfo.YMin =Degree.Second2Degree(rows * dyInSnd); 67 rows = Convert.ToInt32(Math.Round(Degree.Degree2Second(gYMax) / dyInSnd)) + 1; 68 gridInfo.YMax =Degree.Second2Degree( rows * dyInSnd); 69 70 gridInfo.rows = Convert.ToInt32(Math.Round(Degree.Degree2Second(gridInfo.YMax - gridInfo.YMin) / dyInSnd)); 71 gridInfo.cols = Convert.ToInt32(Math.Round(Degree.Degree2Second(gridInfo.XMax - gridInfo.XMin) / dxInSnd)); 72 73 if (envelope.SpatialReference is IProjectedCoordinateSystem) 74 { 75 Geo2Proj(envelope.SpatialReference, gridInfo.XMax, gridInfo.YMax, ref pXMax, ref pYMax); 76 Geo2Proj(envelope.SpatialReference, gridInfo.XMin, gridInfo.YMin, ref pXMin, ref pYMin); 77 gridInfo.XMax = pXMax; gridInfo.XMin = pXMin; 78 gridInfo.YMax = pYMax; gridInfo.YMin = pYMin; 79 } 80 81 return gridInfo; 82 }
(3)图幅命名类
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 public class FrameName 2 { 3 private int scale; 4 /// <summary> 5 /// 6 /// </summary> 7 /// <param name="Scale">分数比例尺的分母(比例尺为1:50000则参数为50000)</param> 8 public FrameName(int Scale) 9 { 10 scale = Scale; 11 } 12 13 public String getFrameName(double longtitude, double latitude) 14 { 15 int baseRowCount=(int)(latitude / 4); 16 char baseChar = Convert.ToChar('A' + baseRowCount); 17 int baseNum = (int)(longtitude / 6) + 31; 18 19 if (scale == 1000000) 20 return String.Format("{0}{1}", baseChar, baseNum); 21 22 char scaleChar = 'E'; //1:50000比例尺的代码为E,其他以后补充 23 switch (scale) 24 { 25 case 500000: 26 scaleChar='B';break; 27 case 250000: 28 scaleChar = 'C';break; 29 case 100000: 30 scaleChar = 'D';break; 31 case 50000: 32 scaleChar='E';break; 33 case 25000: 34 scaleChar = 'F';break; 35 case 10000: 36 scaleChar = 'G';break; 37 case 5000: 38 scaleChar = 'H';break; 39 } 40 41 double dxInSnd, dyInSnd; 42 switch (scale) 43 { 44 case 5000: 45 dxInSnd = 0.03125; 46 dyInSnd = 0.0208333333; 47 break; 48 case 10000: 49 dxInSnd = 0.0625; 50 dyInSnd = 0.0416666667; 51 break; 52 case 25000: 53 dxInSnd = 0.125; 54 dyInSnd = 0.0833333333; 55 break; 56 case 50000: 57 dxInSnd = 0.25; 58 dyInSnd = 0.1666666667; 59 break; 60 case 100000: 61 dxInSnd = 0.5; 62 dyInSnd = 0.3333333333; 63 break; 64 case 250000: 65 dxInSnd = 1.5; 66 dyInSnd = 1; 67 break; 68 case 500000: 69 dxInSnd = 3; 70 dyInSnd = 2; 71 break; 72 default: 73 dxInSnd = 0; 74 dyInSnd = 0; 75 break; 76 } 77 78 if (dxInSnd == dyInSnd && dxInSnd == 0.0) 79 return null; 80 81 int row = (int)(((baseRowCount + 1) * 4 - latitude)/dyInSnd) + 1; 82 //int row = 24-((int)(latitude / dyInSnd)) % 24; 83 int col = (int)((longtitude % 6) /dxInSnd) + 1; 84 85 return String.Format("{0}{1}{2}{3}{4}", baseChar, baseNum, scaleChar, row.ToString().PadLeft(3, '0'), col.ToString().PadLeft(3, '0')); 86 } 87 88 }
(4)调用函数进行分割
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 private void CreateFrame(GridInfo gridInfo, String outDir, int scale) 2 { 3 double xLength = (gridInfo.XMax - gridInfo.XMin) / gridInfo.cols; 4 double yLength = (gridInfo.YMax - gridInfo.YMin) / gridInfo.rows; 5 FrameName frameName = new FrameName(scale); 6 7 for (int i = 0; i < gridInfo.rows; i++) 8 { 9 double yMin = gridInfo.YMin + i * yLength; 10 double yMax = gridInfo.YMin + (i + 1) * yLength; 11 for (int j = 0; j < gridInfo.cols; j++) 12 { 13 double xMin = gridInfo.XMin + j * xLength; 14 double xMax = gridInfo.XMin + (j + 1) * xLength; 15 String pszOutFile; 16 17 if (envelope.SpatialReference is IProjectedCoordinateSystem) 18 { 19 double prjX = 0.0, prjY = 0.0; 20 Proj2Geo(envelope.SpatialReference, (xMax + xMin) / 2, (yMax + yMin) / 2, ref prjX, ref prjY); 21 pszOutFile = frameName.getFrameName(prjX,prjY); 22 } 23 else 24 pszOutFile = frameName.getFrameName((xMax + xMin) / 2, (yMax + yMin) / 2); 25 pszOutFile = outDir + "\" + pszOutFile + ".tif"; 26 CutUpdateImageByAOIGDAL(rasterLyr.FilePath, pszOutFile, xMin, xMax, yMin, yMax, "GTiff"); 27 } 28 } 29 }
通过这4步,一个简易的图幅分割工具就做好了。