在上篇实现了遥感影像的切割，本篇讲切割前的准备。主要分为以下几步：
（1）将影像的投影坐标转换为地理坐标，以便于之后的图幅划分。AE坐标转换函数如下
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）计算包含影像的标准图幅的四至，首先我定义了一个格网类
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 }
根据比例尺计算格网的大小
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）图幅命名类
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）调用函数进行分割
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步，一个简易的图幅分割工具就做好了。
时间： 2024-12-22 17:28:46