MeteoInfo脚本示例：计算每个区域的平均温度

MeteoInfo

最近一个美国人问怎样批量输出GrADS的maskout数据，比如输出每个县的maskout文件，用'Output Map Data'对话框手动来做是很让人抓狂的，因此帮他写了一个脚本来自动处理，为简单起见，行政区域到state一级，也就是把每个state的maskout格点数据输出到一个GrADS格点数据文件中。他想用这些maskout数据来计算每个区域的温度平均等。

需要用到MeteoInfo的最新文件（见置顶帖子）。

其实用MeteoInfo脚本就很容易实现计算每个区域的平均等，因此也顺便写了一个脚本来演示如何实现。用到的美国state一级行政区域shape文件： STATES.SHP (217.13 KB, 下载次数: 62)

2012-11-15 13:58 上传

点击文件名下载附件 ， STATES.SHX (508 Bytes, 下载次数: 41)

2012-11-15 13:58 上传

点击文件名下载附件 ， States.dbf (21.9 KB, 下载次数: 45)

2012-11-15 13:58 上传

点击文件名下载附件 。用到的温度格点示例数据： usgrid.dat (328.56 KB, 下载次数: 65)

2012-11-15 13:59 上传

点击文件名下载附件 。


核心思路是利用GridData的Maskout方法，该方法有几个重载，这里用到Maskout(PolygonShape aPGS)方法，用循环依次用state多边形图元来maskout格点数据，然后再用GridData的Average等方法来计算平均值等。


完整的代码如下：

1. #--------------------------------------------------------        
   
2. # Author: Yaqiang Wang                                          
   
3. # Date: 2012-11-15                                                
   
4. # Purpose: Average the grid data of each polygon
   
5. # Note: Sample                                                   
   
6. #-----------------------------------------------------------
   
7. import clr
   
8. from System.Collections.Generic import *
   
9. clr.AddReference("MeteoInfoC.dll")
   
10. from MeteoInfoC import *
    
11. from MeteoInfoC.Data import *
    
12. from MeteoInfoC.Data.MeteoData import *
    
13. 
    
14. #---- Set directories
    
15. dataDir = "D:\\Temp\"
    
16. 
    
17. #---- Open grid data file - for other data formats: OpenGrADSData, OpenGRIBData ...
    
18. mydata = MeteoDataInfo()
    
19. mydata.OpenSuferGridData(dataDir + "ascii\\usgrid.dat")
    
20. 
    
21. #---- Get states layer
    
22. aLayer = mipy.MapDocument.ActiveMapFrame.MapView.GetLayerFromName("STATES.SHP")
    
23. print aLayer.LayerName
    
24. 
    
25. #---- Mask the grid data by each state polygon and output the average,
    
26. #---- maximum and minimum value of each state
    
27. datagrid = mydata.GetGridData("var")
    
28. i = 0
    
29. for aPolygon in aLayer.ShapeList:
    
30.         maskgrid = datagrid.Maskout(aPolygon)
    
31.         name = aLayer.GetCellValue('STATE_NAME', i)
    
32.         ave = maskgrid.Average()
    
33.         min = maskgrid.GetMinValue()
    
34.         max = maskgrid.GetMaxValue()
    
35.         print name + ', Ave: %.2f, Min: %.2f, Max: %.2f' %(ave, min, max)        
    
36.         i += 1
    
37. 
    
38. print 'Finished!'

复制代码

顺便也贴上批量生成GrADS maskout数据文件的脚本代码：

1. #--------------------------------------------------------        
   
2. # Author: Yaqiang Wang                                          
   
3. # Date: 2012-11-15                                                
   
4. # Purpose: Create multi mask grid data from the polygons  
   
5. # Note: Sample                                                   
   
6. #-----------------------------------------------------------
   
7. import clr
   
8. clr.AddReferenceByPartialName("System")
   
9. clr.AddReferenceByPartialName("System.Drawing")
   
10. from System import *
    
11. from System.Drawing import *
    
12. from System.Collections.Generic import *
    
13. clr.AddReference("MeteoInfoC.dll")
    
14. from MeteoInfoC import *
    
15. from MeteoInfoC.Data import *
    
16. from MeteoInfoC.Data.MeteoData import *
    
17. from MeteoInfoC.Legend import *
    
18. 
    
19. #---- Create a reference grid data
    
20. #---- The parameters are: minimum x, delt x, x number, minimum y, delt y, y number
    
21. xmin = -125.875
    
22. xdelt = 0.25
    
23. xnum = 241
    
24. ymin = 24.125
    
25. ydelt = 0.25
    
26. ynum = 105
    
27. refgrid = GridData(xmin, xdelt, xnum, ymin, ydelt, ynum)
    
28. refgrid.UNDEF = -1
    
29. refgrid.SetValue(1)
    
30. 
    
31. #---- Get states layer
    
32. aLayer = mipy.MapDocument.ActiveMapFrame.MapView.GetLayerFromName("STATES.SHP")
    
33. print aLayer.LayerName
    
34. 
    
35. #---- Set output GrADS data file
    
36. outCtlFile = 'D:\\Temp\\GrADS\\statesmask.ctl'
    
37. outDataFile = 'D:\\Temp\\GrADS\\statesmask.dat'
    
38. 
    
39. #---- Define and set output GrADS data info
    
40. gDataInfo = GrADSDataInfo()
    
41. gDataInfo.FileName = outCtlFile
    
42. gDataInfo.DSET = outDataFile
    
43. gDataInfo.TITLE = 'US states mask out data'
    
44. gDataInfo.DTYPE = 'GRIDDED'
    
45. gDataInfo.XDEF.Type = 'LINEAR'
    
46. gDataInfo.XDEF.XNum = xnum
    
47. gDataInfo.XDEF.XMin = xmin
    
48. gDataInfo.XDEF.XDelt = xdelt
    
49. gDataInfo.YDEF.Type = 'LINEAR'
    
50. gDataInfo.YDEF.YNum = ynum
    
51. gDataInfo.YDEF.YMin = ymin
    
52. gDataInfo.YDEF.YDelt = ydelt
    
53. gDataInfo.ZDEF.Type = 'LINEAR'
    
54. gDataInfo.ZDEF.ZNum = 1
    
55. gDataInfo.ZDEF.SLevel = 1
    
56. gDataInfo.ZDEF.ZDelt = 1
    
57. gDataInfo.TDEF.Type = 'LINEAR'
    
58. gDataInfo.TDEF.TNum = 1
    
59. gDataInfo.TDEF.STime = DateTime(2012, 11, 15)
    
60. gDataInfo.TDEF.TDelt = '1dy'
    
61. i = 0
    
62. for aPolygon in aLayer.ShapeList:
    
63.         name = aLayer.GetCellValue('STATE_NAME', i)
    
64.         aVar = VAR()
    
65.         aVar.VName = 'state' + str(i + 1)
    
66.         aVar.LevelNum = 0
    
67.         aVar.Units = '99'
    
68.         aVar.Description = name
    
69.         gDataInfo.VARDEF.AddVar(aVar)
    
70.         i += 1
    
71. 
    
72. #---- Create mask grid data by each state polygon and write them in GrADS data file
    
73. gDataInfo.CreateDataFile(outDataFile)
    
74. i = 0
    
75. for aPolygon in aLayer.ShapeList:
    
76.         maskgrid = refgrid.Maskout(aPolygon)
    
77.         name = aLayer.GetCellValue('STATE_NAME', i)
    
78.         print str(i) + '  ' + name
    
79.         if i == 0:    #Test one mask grid data by adding and view a layer
    
80.                 maskLayer = mipy.MeteoDataset.Display(maskgrid)
    
81.                 maskLayer.UpdateLegendScheme(LegendType.UniqueValue, 'Data')
    
82.                 maskLayer.LegendScheme.LegendBreaks[0].DrawOutline = False
    
83.                 maskLayer.LegendScheme.LegendBreaks[0].Color = Color.LightGray
    
84.                 maskLayer.LegendScheme.LegendBreaks[1].DrawOutline = False
    
85.                 maskLayer.LegendScheme.LegendBreaks[1].Color = Color.Red
    
86.                 mipy.MapDocument.ActiveMapFrame.MapView.PaintLayers()
    
87.                 mipy.MapDocument.ActiveMapFrame.UpdateLayerNode(maskLayer)
    
88.         
    
89.         gDataInfo.WriteGridData(maskgrid)
    
90.         i += 1
    
91. 
    
92. gDataInfo.CloseDataFile()
    
93. 
    
94. #---- Write GrADS control file
    
95. gDataInfo.WriteGrADSCTLFile()
    
96. 
    
97. print 'Finished!'

复制代码

♂雨已~○

第一个报道学习了，王老师辛苦了！

♂雨已~○

王老师，我用你新更新的文件再执行以前你给的MeteoinfoDemo的时候，下面这几句代码不行了，不知道王老师在最新版的文件中把这个  GridDataPara类弄到哪儿呢？还是给去掉了？用了新更新的文件这个插值功能就不行了！

1.    private void TSMI_StationShaded_Click(object sender, EventArgs e)
   
2.         {
   
3.             this.Cursor = Cursors.WaitCursor;
   
4. 
   
5.             //Read data info
   
6.             MeteoDataInfo aDataInfo = new MeteoDataInfo();
   
7.             string aFile = Application.StartupPath + "\\Sample\\rain_2008072220.csv";
   
8.             aDataInfo.OpenLonLatData(aFile);
   
9. 
   
10.             //Get station data
    
11.             StationData stationData = aDataInfo.GetStationData("Precipitation6h");
    
12. 
    
13.             //Interpolate
    
14.             GridDataPara aGDP = new GridDataPara();
    
15.             aGDP.dataExtent.minX = 60;
    
16.             aGDP.dataExtent.maxX = 140;
    
17.             aGDP.dataExtent.minY = -20;
    
18.             aGDP.dataExtent.maxY = 60;
    
19.             aGDP.xNum = 80;
    
20.             aGDP.yNum = 80;
    
21.             GridInterpolation gridInterp = new GridInterpolation();
    
22.             gridInterp.GridDataParaV = aGDP;
    
23. 
    
24.             gridInterp.GridInterMethodV = GridInterMethod.IDW_Radius;
    
25.             gridInterp.Radius = 2;
    
26.             gridInterp.MinPointNum = 1;
    
27. 
    
28.             double[] X = new double[1];
    
29.             double[] Y = new double[1];
    
30.             ContourDraw.CreateGridXY(gridInterp.GridDataParaV, ref X, ref Y);
    
31.             double[,] S = stationData.Data;
    
32.             S = ContourDraw.FilterDiscreteData_Radius(S, gridInterp.Radius,
    
33.                 gridInterp.GridDataParaV.dataExtent, stationData.UNDEF);
    
34.             GridData gridData = ContourDraw.InterpolateDiscreteData_Radius(S,
    
35.                 X, Y, gridInterp.MinPointNum, gridInterp.Radius, stationData.UNDEF);
    
36. 
    
37.             //Create legend scheme
    
38.             bool hasNoData = true;
    
39.             LegendScheme aLS = LegendManage.CreateLegendSchemeFromGridData(gridData, LegendType.GraduatedColor,
    
40.                 ShapeTypes.Polygon, ref hasNoData);
    
41.             ((PolygonBreak)aLS.breakList[0]).DrawFill = false;
    
42. 
    
43.             //Create layer
    
44.             VectorLayer aLayer = new VectorLayer(ShapeTypes.Polygon);
    
45.             aLayer = DrawMeteoData.CreateShadedLayer(gridData, aLS, "Rain", "Rain");
    
46.             aLayer.IsMaskout = true;
    
47. 
    
48.             //Add layer
    
49.             layersLegend1.ActiveMapFrame.AddLayer(aLayer);
    
50.             layersLegend1.ActiveMapFrame.MoveLayer(aLayer.Handle, 0);
    
51.             layersLegend1.Refresh();
    
52. 
    
53.             //Change title of the layout
    
54.             LayoutGraphic aTitle = mapLayout1.GetTexts()[0];
    
55.             aTitle.SetLabelText("MeteoInfo Class Library Demo - Station Shaded Layer");
    
56. 
    
57.             //Add or change the legend in layout
    
58.             LayoutLegend aLegend;
    
59.             if (mapLayout1.GetLegends().Count > 0)
    
60.                 aLegend = mapLayout1.GetLegends()[0];
    
61.             else
    
62.                 aLegend = mapLayout1.AddLegend(650, 100);
    
63.             aLegend.LegendStyle = LegendStyleEnum.Bar_Vertical;
    
64.             aLegend.LegendLayer = aLayer;
    
65.             if (tabControl1.SelectedIndex == 1)
    
66.                 mapLayout1.PaintGraphics();
    
67. 
    
68.             this.Cursor = Cursors.Default;

复制代码

♂雨已~○

王老师，我尝试着用你这种方法来先将数据maskout后作图，但是出现了只有左边图例，没有图形，这是为什么？下面是我的代码，请王老师帮忙看看~~~~

1.            //Create a legend scheme
   
2.            MeteoDataInfo aDataInfo = new MeteoDataInfo();
   
3. 
   
4.            //Open GrADS data file
   
5.            string aFile = Application.StartupPath + "\\Sample\\model.ctl";
   
6.            aDataInfo.OpenGrADSData(aFile);
   
7. 
   
8.            //Get GridData
   
9.            GridData press = aDataInfo.GetGridData("PS");
   
10. 
    
11.           aFile = Application.StartupPath + "\\Map\\四川重庆.shp";
    
12.           MapLayer maskLayer = MapDataManage.OpenLayer(aFile);
    
13.           layersLegend1.ActiveMapFrame.AddLayer(maskLayer);
    
14.           layersLegend1.Refresh();
    
15.               
    
16. 
    
17.            //Create a legend scheme
    
18.             
    
19.            GridData maskpress = press.Maskout(mapView1,"四川重庆.shp");
    
20. 
    
21.            bool hasUndefData = false;
    
22.            LegendScheme aLS = LegendManage.CreateLegendSchemeFromGridData(maskpress,
    
23.                        LegendType.GraduatedColor, ShapeTypes.Polygon, ref hasUndefData);
    
24. 
    
25.            //Create a contour layer
    
26.            VectorLayer aLayer = DrawMeteoData.CreateShadedLayer(maskpress, aLS, "Shaded_PS", "PS");
    
27.            
    
28.            //Add layer
    
29.            layersLegend1.ActiveMapFrame.AddLayer(aLayer);
    
30.            layersLegend1.ActiveMapFrame.MoveLayer(aLayer.Handle, 0);
    
31.            layersLegend1.Refresh();
    
32. 
    
33.            //Change title of the layout
    
34.            LayoutGraphic aTitle = mapLayout1.GetTexts()[0];
    
35.            aTitle.SetLabelText("MeteoInfo Class Library Demo - Shaded Layer");
    
36. 
    
37.            //Add or change the legend in layout
    
38.            LayoutLegend aLegend;
    
39.            if (mapLayout1.GetLegends().Count > 0)
    
40.                aLegend = mapLayout1.GetLegends()[0];
    
41.            else
    
42.                aLegend = mapLayout1.AddLegend(650, 100);
    
43.            aLegend.LegendStyle = LegendStyleEnum.Bar_Vertical;
    
44.            aLegend.LegendLayer = aLayer;
    
45.            if (tabControl1.SelectedIndex == 1)
    
46.                mapLayout1.PaintGraphics();

复制代码

MeteoInfo

♂雨已~○ 发表于 2012-11-15 19:18

                               
登录/注册后可看大图

王老师，我用你新更新的文件再执行以前你给的MeteoinfoDemo的时候，下面这几句代码不行了，不知道王老师在最 ...

改了一些类的名称，你在http://www.meteothinker.com/上下载最新的Demo程序看看相关代码就知道了。

MeteoInfo

♂雨已~○ 发表于 2012-11-15 19:30

                               
登录/注册后可看大图

王老师，我尝试着用你这种方法来先将数据maskout后作图，但是出现了只有左边图例，没有图形，这是为什么？下 ...

model.ctl数据分辨率很粗，在重庆市区域也就一个点，当然做不出等值线图来了。

♂雨已~○

MeteoInfo 发表于 2012-11-15 20:49

                               
登录/注册后可看大图

model.ctl数据分辨率很粗，在重庆市区域也就一个点，当然做不出等值线图来了。

奥，明白了，谢谢王老师，我用站点资料就行了~~~

denghaijun

学习了，谢谢王老师！

卉卉

MeteoInfo 发表于 2012-11-15 20:46

                               
登录/注册后可看大图

改了一些类的名称，你在http://www.meteothinker.com/上下载最新的Demo程序看看相关代码就知道了。

王老师，要是初学meteinfo的脚本应该从什么地方开始？

MeteoInfo

卉卉 发表于 2012-12-13 20:40

                               
登录/注册后可看大图

王老师，要是初学meteinfo的脚本应该从什么地方开始？

学习IronPython， MeteoInfo类库，论坛里的脚本帖子