• AE + GDAL实现影像按标准图幅分割（下）

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

    　　（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         }
    View Code

    　　（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     }
    View Code

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

     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         }
    View Code

    　　（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     }
    View Code

    　　（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         }
    View Code

    　　通过这4步，一个简易的图幅分割工具就做好了。
  • 相关阅读:
    Shell-01-脚本开头实现自动添加注释
    Linux中通过SHELL发送邮件
    linux服务器修改密码登录Failed to restart ssh.service: Unit ssh.service not found
    ffmpeg+java实现五秒钟剪辑80个视频
    Vue学习-watch 监听用法
    springboot添加定时任务
    Spring异常：java.lang.NoClassDefFoundError: org/springframework/core/OrderComparator$OrderSourceProvider
    多线程实战-龟兔赛跑
    Git分支管理（二）
    android studio bug : aidl is missing 解决方案
  • 原文地址：https://www.cnblogs.com/sunnyeveryday/p/4518389.html
Copyright © 2020-2023  润新知