地统计分析在气象领域的应用

由于工作需要，要求使用AE实现对某些气象观测要素（如气温、雨量）等进行IDW插值，经过这段时间的努力，基本功能已经实现。在此感谢一些网上的技术牛人，谢谢他们无私的分享（搜索是件快乐的事情），同时也要感谢自己付出的努力（智慧、查找资料、耐心）。实现过程大致如下：

IDW插值
剪裁（单波段影像）
颜色渲染
出透明图

下面记录了实现的主要代码，毕竟好记性不如烂笔头。

（一）IDW空间插值的实现

 1 private ESRI.ArcGIS.Geodatabase.IGeoDataset CreateRaster(string filedName, string pPath)
 2         {
 3             IMap pMap = this._mapControl.Map;
 4             IInterpolationOp3 pInterpolationOp = new RasterInterpolationOpClass(); //IInterpolationOp目前已被IInterpolationOp3所取代
 5             //定义工作空间
 6             IWorkspaceFactory pWorkspaceFactory = new ShapefileWorkspaceFactory();
 7             string pFolder = System.IO.Path.GetDirectoryName(pPath);
 8             string pFileName = System.IO.Path.GetFileName(pPath);
 9             IWorkspace pWorkspace = pWorkspaceFactory.OpenFromFile(pFolder, 0);  //设定地理处理的工作空间
10             IFeatureWorkspace pFeatureWorkspace = pWorkspace as IFeatureWorkspace;
11             IFeatureClass oFeatureClass = pFeatureWorkspace.OpenFeatureClass(pFileName);
12             IFeatureClassDescriptor pFCDescriptor = new FeatureClassDescriptorClass();
13             pFCDescriptor.Create(oFeatureClass, null, filedName);
14             IRasterRadius pRadius = new RasterRadiusClass();     //设置搜索半径
15 
16             object objectMaxDistance = null;
17             object objectbarrier = null;
18             object missing = Type.Missing;
19             pRadius.SetVariable(12, ref missing);  //这里设置不同的值，出图的效果（如颜色的深浅等）就不同
20 
21             //ESRI.ArcGIS.Geodatabase.IWorkspace pShpWorkspace = pWorkspaceFactory.OpenFromFile(tempPath, 0);
22 
23             object dCellSize = 0.013;   //设置栅格精度（输出栅格像元大小）
24             object snapRasterData = Type.Missing;
25             IRasterAnalysisEnvironment pEnv = pInterpolationOp as IRasterAnalysisEnvironment;
26             pEnv.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, ref dCellSize);
27             pEnv.OutSpatialReference = pMap.SpatialReference;
28             //执行插置，产生结果（栅格数据集）
29             IGeoDataset poutGeoDataset = pInterpolationOp.IDW((IGeoDataset)pFCDescriptor, 2, pRadius, ref missing);  //这里参数2为默认的权重值
30             //控制周围点对于内插值的重要性的距离指数。幂值越高，远数据点的影响会越小。
31             IRaster pOutRaster = poutGeoDataset as IRaster;
32 
33             IRasterLayer pOutRasLayer = new RasterLayer();
34             pOutRasLayer.CreateFromRaster(pOutRaster);  //生成栅格数据图层
35             this._mapControl.AddLayer(pOutRasLayer);   //添加生成的图层
36 
37             return poutGeoDataset;  //返回一个GeoDataSet
38             
39         }

（二）shp图层剪裁栅格数据

 1 private IRaster ShpLayerClipRaster(IFeatureLayer featureLayer, IRaster raster)
 2         {
 3             IFeatureLayer pFeaLyr;
 4             IRasterAnalysisEnvironment pEnv;
 5             IGeoDataset pTempDS;
 6             pFeaLyr = featureLayer;
 7             pTempDS = pFeaLyr.FeatureClass as IGeoDataset;
 8             IConversionOp pConOp = new RasterConversionOpClass();
 9             pEnv = pConOp as IRasterAnalysisEnvironment;
10             IRasterProps pProp;
11             pProp = raster as IRasterProps;
12             object cellSize = 0.01;
13             pEnv.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, ref cellSize);
14             IRasterConvertHelper rsh = new RasterConvertHelperClass();
15             IRaster tempRaster;
16             tempRaster = rsh.ToRaster1(pTempDS, "Grid", pEnv);
17 
18             IRaster pOutRaster;
19             IExtractionOp pExtrOp = new RasterExtractionOpClass();
20             pOutRaster = (IRaster)pExtrOp.Raster((IGeoDataset)raster, (IGeoDataset)tempRaster);
21 
22             //pTempDS = null;
23             //pConOp = null;
24             //pFeaLyr = null;   
25 
26             return pOutRaster;
27         }

（三）栅格数据颜色分级渲染

实现思路如下：

1.定义渲染的一系列接口

2.判断图像是否建立了直方图，如果没有则进行创建

3.定义颜色序列，为渲染提供渲染的方案（自定义色带颜色包括透明效果色以及范围值区间）

4.调用Render方法进行渲染

  1 private void classifyColorToRaster(IRasterLayer pRasterLayer, int classCount, int size)
  2         {
  3             IRasterClassifyColorRampRenderer classifyRender = new RasterClassifyColorRampRenderer();
  4             IRasterRenderer rasterRender = classifyRender as IRasterRenderer;
  5 
  6             IRaster pRaster = pRasterLayer.Raster;
  7             IRasterBandCollection pRBandCol = pRaster as IRasterBandCollection;
  8             IRasterBand pRBand = pRBandCol.Item(0);
  9 
 10             if (pRBand.Histogram == null)   //判断栅格是否计算直方图，没有的话要进行计算
 11             {
 12                 pRBand.ComputeStatsAndHist();
 13             }
 14             rasterRender.Raster = pRaster;
 15             classifyRender.ClassCount = classCount;
 16             
 17             for (int i = 0; i < classCount; i++)
 18             {
 19                 if (i == 0)
 20                 {
 21                     classifyRender.set_Break(i, double.Parse("1"));
 22                 }
 23                 if (i == 1)
 24                 {
 25                     classifyRender.set_Break(i, double.Parse("251"));
 26                 }
 27                 if (i == 2)
 28                 {
 29                     classifyRender.set_Break(i, double.Parse("501"));
 30                 }
 31                 if (i == 3)
 32                 {
 33                     classifyRender.set_Break(i, double.Parse("751"));
 34                 }
 35                 if (i == 4)
 36                 {
 37                     classifyRender.set_Break(i, double.Parse("1001"));
 38                 }
 39                 if (i == 5)
 40                 {
 41                     classifyRender.set_Break(i, double.Parse("1251"));
 42                 }
 43                 if (i == 6)
 44                 {
 45                     classifyRender.set_Break(i, double.Parse("1501"));
 46                 }
 47                 if (i == 7)
 48                 {
 49                     classifyRender.set_Break(i, double.Parse("1751"));
 50                 }
 51                 if (i == 8)
 52                 {
 53                     classifyRender.set_Break(i, double.Parse("2001"));
 54                 }
 55                 if (i == 9)
 56                 {
 57                     classifyRender.set_Break(i, double.Parse("2251"));
 58                 }
 59             }
 60             
 61             rasterRender.Update();
 62             /*
 63             IRgbColor pFromColor = new RgbColor() as IRgbColor;
 64             pFromColor.Red = formColorRed;
 65             pFromColor.Green = formColorGreen;
 66             pFromColor.Blue = formColorBlue;
 67             IRgbColor pToColor = new RgbColor() as IRgbColor;
 68             pToColor.Red = toColorRed;
 69             pToColor.Green = toColorGreen;
 70             pToColor.Blue = toColorBlue;
 71             */
 72             IAlgorithmicColorRamp colorRamp = new AlgorithmicColorRamp() as IAlgorithmicColorRamp;
 73             colorRamp.Size = size;
 74             //colorRamp.FromColor = pFromColor;
 75             //colorRamp.ToColor = pToColor;
 76             bool createColorRamp;
 77 
 78             colorRamp.CreateRamp(out createColorRamp);
 79 
 80             IFillSymbol fillSymbol = new SimpleFillSymbol() as IFillSymbol;
 81             for (int i = 0; i < classifyRender.ClassCount; i++)
 82             {
 83                 if(i == 0)
 84                 {
 85                     ESRI.ArcGIS.Display.IRgbColor tc = new ESRI.ArcGIS.Display.RgbColor();
 86                     tc.Red = 255;
 87                     tc.Green = 255;
 88                     tc.Blue = 128;
 89                     fillSymbol.Color = tc as IColor;
 90                     classifyRender.set_Symbol(i, fillSymbol as ISymbol);
 91                     classifyRender.set_Label(i, classifyRender.get_Break(i).ToString("0") + "-" + "250");  //这里设置0是要整数显示
 92                 }
 93                 if (i == 1)
 94                 {
 95                     ESRI.ArcGIS.Display.IRgbColor tc = new ESRI.ArcGIS.Display.RgbColor();
 96                     tc.Red = 128;
 97                     tc.Green = 255;
 98                     tc.Blue = 128;
 99                     fillSymbol.Color = tc as IColor;
100                     classifyRender.set_Symbol(i, fillSymbol as ISymbol);
101                     classifyRender.set_Label(i, classifyRender.get_Break(i).ToString("0") + "-" + "500");  //这里设置0是要整数显示
102                 }
103                 if (i == 2)
104                 {
105                     ESRI.ArcGIS.Display.IRgbColor tc = new ESRI.ArcGIS.Display.RgbColor();
106                     tc.Red = 128;
107                     tc.Green = 255;
108                     tc.Blue = 255;
109                     fillSymbol.Color = tc as IColor;
110                     classifyRender.set_Symbol(i, fillSymbol as ISymbol);
111                     classifyRender.set_Label(i, classifyRender.get_Break(i).ToString("0") + "-" + "750");  //这里设置0是要整数显示
112                 }
113                 if (i == 3)
114                 {
115                     ESRI.ArcGIS.Display.IRgbColor tc = new ESRI.ArcGIS.Display.RgbColor();
116                     tc.Red = 0;
117                     tc.Green = 128;
118                     tc.Blue = 255;
119                     fillSymbol.Color = tc as IColor;
120                     classifyRender.set_Symbol(i, fillSymbol as ISymbol);
121                     classifyRender.set_Label(i, classifyRender.get_Break(i).ToString("0") + "-" + "1000");  //这里设置0是要整数显示
122                 }
123                 if (i == 4)
124                 {
125                     ESRI.ArcGIS.Display.IRgbColor tc = new ESRI.ArcGIS.Display.RgbColor();
126                     tc.Red = 255;
127                     tc.Green = 128;
128                     tc.Blue = 192;
129                     fillSymbol.Color = tc as IColor;
130                     classifyRender.set_Symbol(i, fillSymbol as ISymbol);
131                     classifyRender.set_Label(i, classifyRender.get_Break(i).ToString("0") + "-" + "1250");  //这里设置0是要整数显示
132                 }
133                 if (i == 5)
134                 {
135                     ESRI.ArcGIS.Display.IRgbColor tc = new ESRI.ArcGIS.Display.RgbColor();
136                     tc.Red = 255;
137                     tc.Green = 0;
138                     tc.Blue = 0;
139                     fillSymbol.Color = tc as IColor;
140                     classifyRender.set_Symbol(i, fillSymbol as ISymbol);
141                     classifyRender.set_Label(i, classifyRender.get_Break(i).ToString("0") + "-" + "1500");  //这里设置0是要整数显示
142                 }
143                 if (i == 6)
144                 {
145                     ESRI.ArcGIS.Display.IRgbColor tc = new ESRI.ArcGIS.Display.RgbColor();
146                     tc.Red = 255;
147                     tc.Green = 0;
148                     tc.Blue = 255;
149                     fillSymbol.Color = tc as IColor;
150                     classifyRender.set_Symbol(i, fillSymbol as ISymbol);
151                     classifyRender.set_Label(i, classifyRender.get_Break(i).ToString("0") + "-" + "1750");  //这里设置0是要整数显示
152                 }
153                 if (i == 7)
154                 {
155                     ESRI.ArcGIS.Display.IRgbColor tc = new ESRI.ArcGIS.Display.RgbColor();
156                     tc.Red = 255;
157                     tc.Green = 128;
158                     tc.Blue = 64;
159                     fillSymbol.Color = tc as IColor;
160                     classifyRender.set_Symbol(i, fillSymbol as ISymbol);
161                     classifyRender.set_Label(i, classifyRender.get_Break(i).ToString("0") + "-" + "2000");  //这里设置0是要整数显示
162                 }
163                 if (i == 8)
164                 {
165                     ESRI.ArcGIS.Display.IRgbColor tc = new ESRI.ArcGIS.Display.RgbColor();
166                     tc.Red = 0;
167                     tc.Green = 64;
168                     tc.Blue = 128;
169                     fillSymbol.Color = tc as IColor;
170                     classifyRender.set_Symbol(i, fillSymbol as ISymbol);
171                     classifyRender.set_Label(i, classifyRender.get_Break(i).ToString("0") + "-" + "2250");  //这里设置0是要整数显示
172                 }
173                 if (i == 9)
174                 {
175                     ESRI.ArcGIS.Display.IRgbColor tc = new ESRI.ArcGIS.Display.RgbColor();
176                     tc.Red = 255;
177                     tc.Green = 255;
178                     tc.Blue = 0;
179                     fillSymbol.Color = tc as IColor;
180                     classifyRender.set_Symbol(i, fillSymbol as ISymbol);
181                     classifyRender.set_Label(i, classifyRender.get_Break(i).ToString("0") + "-" + "2500");  //这里设置0是要整数显示
182                 }
183             }
184             pRasterLayer.Renderer = rasterRender;
185             //*********************************使栅格图层呈现透明效果*******************************************
186             //ILayerEffects pLayerEffects = pRasterLayer as ILayerEffects;    
187             //pLayerEffects.Transparency = 50;
188             //*************************************************************************************************
189 
190 
191             this._mapControl.AddLayer(pRasterLayer);
192             this._mapControl.ActiveView.Refresh();
193         }

（四）导出透明效果图片（PNG格式）

 1 private void ExportMapByAbsPath(IActiveView activeView, string path)     //导出透明效果图片
 2         {
 3             IExport export = new ExportPNGClass();
 4             IExportImage exportImg;
 5             IExportPNG exportPNG;
 6             int resolution = 96;
 7             export.Resolution = resolution;
 8             exportImg = export as IExportImage;
 9             exportImg.ImageType = esriExportImageType.esriExportImageTypeTrueColor;  //定义输出图片颜色模式为24位真彩色
10             exportPNG = export as IExportPNG;
11 
12             //exportImg.BackgroundColor.RGB = exportPNG.TransparentColor.RGB;
13             ESRI.ArcGIS.Display.IRgbColor tc = new ESRI.ArcGIS.Display.RgbColor();
14             tc.Red = 255;
15             tc.Green = 255;
16             tc.Blue = 0;
17             exportImg.BackgroundColor = tc as IColor;   //设置输出图片的背景色
18             exportPNG.TransparentColor = tc as IColor;  //设置输出图片的透明色
19             tagRECT exportRect = activeView.ExportFrame;
20             IEnvelope envelop = new EnvelopeClass();
21             envelop.PutCoords(exportRect.left, exportRect.top, exportRect.right, exportRect.bottom);
22             export.PixelBounds = envelop;
23             export.ExportFileName = path;
24             int hDC = export.StartExporting();
25             IEnvelope pVisbounds = null;
26             ITrackCancel ptrac = null;
27             //导出
28             activeView.Output(hDC, (int)export.Resolution, ref exportRect, pVisbounds, ptrac);
29             //结束导出
30             export.FinishExporting();
31             //清理导出类
32             export.Cleanup();
33         }

（五）调用上述方法

 1 private void exportMapBtn_Click(object sender, EventArgs e)    //实现了shp文件剪裁插值结果（栅格）
 2         {
 3             string filedName = "RES1_4M0_I";
 4             string pPath = System.Windows.Forms.Application.StartupPath + @"china_basic_map省会城市.shp";
 5             CreateRaster(filedName, pPath);
 6             IRasterLayer pRasterLayer = this._mapControl.get_Layer(0) as RasterLayer;
 7             IFeatureLayer featureLayer = this._mapControl.get_Layer(2) as IFeatureLayer;
 8             IRaster raster = pRasterLayer.Raster;
 9             IRaster resultRaster = ShpLayerClipRaster(featureLayer, raster);
10             IRasterLayer resultRasterLayer = new RasterLayer();
11             resultRasterLayer.CreateFromRaster(resultRaster);
12             classifyColorToRaster(resultRasterLayer, 10, 10);
13             //***************************该删除方式有待改善*************************************
14             this._mapControl.DeleteLayer(1);
15             this._mapControl.DeleteLayer(1);
16             this._mapControl.DeleteLayer(1);
17             //*********************************************************************************
18             ExportMapByAbsPath(this._mapControl.ActiveView, "d:\abc.png");  //输出png格式图片
19             this._mapControl.DeleteLayer(0);
20         }

（五）目前任然存在的问题以及下一步的需求，后续有待解决

   1.License偶尔报错（概率很小）
   2.自定义色带时，初始值设置为0时，颜色渲染有问题，只要比0大便不会出现此问题
   3.剪裁时需先调整范围，才能使剪裁范围最全面，效果最优
   4.清理缓存，提高运行效率
   5.Janus Form为试用版，最好要找到破解版
   6.代码优化（组织思维、重构）
   7.数据接入方式
   8.定时出图（分钟级、小时级、日级、周级、月级、旬级、年级）