地統計分析在氣象領域的應用

因爲工做須要，要求使用AE實現對某些氣象觀測要素（如氣溫、雨量）等進行IDW插值，通過這段時間的努力，基本功能已經實現。在此感謝一些網上的技術牛人，謝謝他們無私的分享（搜索是件快樂的事情），同時也要感謝本身付出的努力（智慧、查找資料、耐心）。實現過程大體以下：

1. IDW插值
2. 剪裁（單波段影像）
3. 顏色渲染
4. 出透明圖

下面記錄了實現的主要代碼，畢竟好記性不如爛筆頭。

（一）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.定時出圖（分鐘級、小時級、日級、周級、月級、旬級、年級）