MeteoInfo脚本示例：GrADS to netCDF

MeteoInfo

这里给出一个将GrADS数据文件转为netCDF数据文件的脚本示例程序，其它格式数据转netCDF可以参考：

1. #-----------------------------------------------------
   
2. # Author: Yaqiang Wang
   
3. # Date: 2015-3-12
   
4. # Purpose: Convert CUACE grads data to netCDF (CUACE/Dust)
   
5. # Note: Sample
   
6. #-----------------------------------------------------
   
7. print 'Loading classes...'
   
8. from org.meteoinfo.data import GridData
   
9. from org.meteoinfo.data import DataMath
   
10. from org.meteoinfo.data.meteodata import MeteoDataInfo
    
11. from org.meteoinfo.geoprocess.analysis import ResampleMethods
    
12. from org.meteoinfo.data.meteodata.netcdf import NetCDFDataInfo
    
13. from org.meteoinfo.projection import ProjectionInfo
    
14. from org.meteoinfo.projection import ProjectionNames
    
15. from org.meteoinfo.projection import KnownCoordinateSystems
    
16. from ucar.nc2 import NetcdfFileWriter
    
17. from ucar.nc2 import Attribute
    
18. from ucar.ma2 import DataType
    
19. from ucar.ma2 import Array
    
20. import os.path
    
21. import jarray
    
22. import datetime
    
23. from java.util import Date
    
24. from java.text import SimpleDateFormat
    
25. 
    
26. #Set date
    
27. year = 2014
    
28. month = 4
    
29. day = 23
    
30. hour = 0
    
31. sdate = datetime.datetime(year, month, day, hour)
    
32. print sdate
    
33. #Set directory
    
34. dataDir = 'U:/data/cuace_dust/dust_example/2014_case'
    
35. outDir = dataDir
    
36. infn = os.path.join(dataDir, 'dust_post_'+ sdate.strftime('%Y%m%d%H') + '.ctl')
    
37. outfn = os.path.join(dataDir, 'WMO_SDS-WAS_Asian_Center_Model_Forecasting_CUACE-Dust_CMA_' \
    
38.     + sdate.strftime('%Y-%m-%d') + '.nc')
    
39. #Set output X/Y coordinates and projection
    
40. toProjInfo = KnownCoordinateSystems.geographic.world.WGS1984
    
41. sx = 70.0
    
42. xnum = 161
    
43. sy = 20
    
44. ynum = 71
    
45. delta = 0.5
    
46. xlist = []
    
47. ylist = []
    
48. for i in range(0, xnum):
    
49.     xlist.append(sx + delta * i)
    
50. for i in range(0, ynum):
    
51.     ylist.append(sy + delta * i)
    
52. X = jarray.array(xlist, 'd')
    
53. Y = jarray.array(ylist, 'd')
    
54. 
    
55. #Read GrADS data file
    
56. print 'Open GrADS data file...'
    
57. mdi = MeteoDataInfo()
    
58. mdi.openGrADSData(infn)
    
59. dataInfo = mdi.getDataInfo()
    
60. dataInfo.setBigEndian(True)
    
61. fromProjInfo = mdi.getProjectionInfo()
    
62. tnum = dataInfo.getTimeNum()
    
63. mvalue = dataInfo.getMissingValue()
    
64. 
    
65. #Set output nc data file
    
66. print 'Create output NC file: ' + outfn
    
67. ncfile = NetcdfFileWriter.createNew(NetcdfFileWriter.Version.netcdf3, outfn)
    
68. #Add dimensions
    
69. print 'Add dimensions...'
    
70. xDim = ncfile.addDimension(None, 'lon', xnum)
    
71. yDim = ncfile.addDimension(None, 'lat', ynum)
    
72. tDim = ncfile.addDimension(None, 'time', tnum)
    
73. 
    
74. #Add global attributes
    
75. print 'Add global attributes...'
    
76. ncfile.addGroupAttribute(None, Attribute('Conventions', 'CF-1.6'))
    
77. ncfile.addGroupAttribute(None, Attribute('Title', 'Sand and dust storm forecasting'))
    
78. ncfile.addGroupAttribute(None, Attribute('Model', 'CUACE/Dust'))
    
79. ncfile.addGroupAttribute(None, Attribute('Center', 'WMO SDS-WAS Asian Center'))
    
80. ncfile.addGroupAttribute(None, Attribute('Agency', 'China Meteorological Administration'))
    
81. 
    
82. #Add variables
    
83. xvar = ncfile.addVariable(None, 'lon', DataType.FLOAT, [xDim])
    
84. xvar.addAttribute(Attribute('units', 'degrees_east'))
    
85. xvar.addAttribute(Attribute('long_name', 'Longitude'))
    
86. xvar.addAttribute(Attribute('standard_name', 'longitude'))
    
87. xvar.addAttribute(Attribute('axis', 'X'))
    
88. yvar = ncfile.addVariable(None, 'lat', DataType.FLOAT, [yDim])
    
89. yvar.addAttribute(Attribute('units', 'degrees_north'))
    
90. yvar.addAttribute(Attribute('long_name', 'Latitude'))
    
91. yvar.addAttribute(Attribute('standard_name', 'latitude'))
    
92. yvar.addAttribute(Attribute('axis', 'Y'))
    
93. tvar = ncfile.addVariable(None, 'time', DataType.INT, [tDim])
    
94. tvar.addAttribute(Attribute('units', 'hours since 1900-01-01 00:00:0.0'))
    
95. tvar.addAttribute(Attribute('long_name', 'Time'))
    
96. tvar.addAttribute(Attribute('standart_name', 'time'))
    
97. tvar.addAttribute(Attribute('axis', 'T'))
    
98. #Data variables
    
99. vnames = ['load','con','dry','wet','aod']
    
100. varlist = []
     
101. var = ncfile.addVariable(None, 'LOAD_DUST', DataType.FLOAT, [tDim, yDim, xDim])
     
102. var.addAttribute(Attribute('long_name', 'Dust load'))
     
103. var.addAttribute(Attribute('units', 'kg/m2'))
     
104. var.addAttribute(Attribute('missing_value', -9999.0))
     
105. varlist.append(var)
     
106. var = ncfile.addVariable(None, 'SCONC_DUST', DataType.FLOAT, [tDim, yDim, xDim])
     
107. var.addAttribute(Attribute('long_name', 'Surface dust concentration'))
     
108. var.addAttribute(Attribute('units', 'ug/m3'))
     
109. var.addAttribute(Attribute('missing_value', -9999.0))
     
110. varlist.append(var)
     
111. var = ncfile.addVariable(None, 'DDEPO_DUST', DataType.FLOAT, [tDim, yDim, xDim])
     
112. var.addAttribute(Attribute('long_name', '3-hour accumulated dry deposition'))
     
113. var.addAttribute(Attribute('units', 'kg/m2'))
     
114. var.addAttribute(Attribute('missing_value', -9999.0))
     
115. varlist.append(var)
     
116. var = ncfile.addVariable(None, 'WDEPO_DUST', DataType.FLOAT, [tDim, yDim, xDim])
     
117. var.addAttribute(Attribute('long_name', '3-hour accumulated wet deposition'))
     
118. var.addAttribute(Attribute('units', 'kg/m2'))
     
119. var.addAttribute(Attribute('missing_value', -9999.0))
     
120. varlist.append(var)
     
121. var = ncfile.addVariable(None, 'AOD550_DUST', DataType.FLOAT, [tDim, yDim, xDim])
     
122. var.addAttribute(Attribute('long_name', 'Dust optical depth at 550nm'))
     
123. var.addAttribute(Attribute('units', '-'))
     
124. var.addAttribute(Attribute('missing_value', -9999.0))
     
125. varlist.append(var)
     
126. 
     
127. #Write nc file
     
128. ncfile.create()
     
129. #Write x,y,z,t variables
     
130. print 'Write x variable...'
     
131. shape = jarray.array([xnum], 'i')
     
132. data = Array.factory(DataType.FLOAT, shape)
     
133. for i in range(0,xnum):
     
134.     data.set(i, X)
     
135. ncfile.write(xvar, data)
     
136. 
     
137. print 'Write y variable...'
     
138. shape = jarray.array([ynum], 'i')
     
139. data = Array.factory(DataType.FLOAT, shape)
     
140. for i in range(0,ynum):
     
141.     data.set(i, Y)
     
142. ncfile.write(yvar, data)
     
143. 
     
144. print 'Write time variable...'
     
145. format = SimpleDateFormat('yyyy-MM-dd')
     
146. sdate = format.parse('1900-01-01')
     
147. tvalues = dataInfo.getTimeValues(sdate, 'hours')
     
148. shape = jarray.array([tnum], 'i')
     
149. data = Array.factory(DataType.INT, shape)
     
150. for i in range(0,tnum):
     
151.     data.set(i, tvalues)
     
152. ncfile.write(tvar, data)
     
153. 
     
154. #Write data variables
     
155. print 'Write data variable...'
     
156. for vname, var in zip(vnames, varlist):
     
157.     for t in range(0, tnum):
     
158.         print 'Time: ' + str(t + 1)
     
159.         mdi.setTimeIndex(t)
     
160.         gData = mdi.getGridData(vname)
     
161.         ngData = gData.project(fromProjInfo, toProjInfo, X, Y, ResampleMethods.Bilinear)
     
162.         origin = jarray.array([t, 0, 0], 'i')
     
163.         ncfile.write(var, origin, NetCDFDataInfo.gridToArray3D(ngData))
     
164. 
     
165. #Close nc file
     
166. ncfile.flush()
     
167. ncfile.close()
     
168. 
     
169. print 'Finished'

#-----------------------------------------------------<br /> # Author: Yaqiang Wang<br /> # Date: 2015-3-12<br /> # Purpose: Convert CUACE grads data to netCDF (CUACE/Dust)<br /> # Note: Sample<br /> #-----------------------------------------------------<br /> print 'Loading classes...'<br /> from org.meteoinfo.data import GridData<br /> from org.meteoinfo.data import DataMath<br /> from org.meteoinfo.data.meteodata import MeteoDataInfo<br /> from org.meteoinfo.geoprocess.analysis import ResampleMethods<br /> from org.meteoinfo.data.meteodata.netcdf import NetCDFDataInfo<br /> from org.meteoinfo.projection import ProjectionInfo<br /> from org.meteoinfo.projection import ProjectionNames<br /> from org.meteoinfo.projection import KnownCoordinateSystems<br /> from ucar.nc2 import NetcdfFileWriter<br /> from ucar.nc2 import Attribute<br /> from ucar.ma2 import DataType<br /> from ucar.ma2 import Array<br /> import os.path<br /> import jarray<br /> import datetime<br /> from java.util import Date<br /> from java.text import SimpleDateFormat<br /> <br /> #Set date<br /> year = 2014<br /> month = 4<br /> day = 23<br /> hour = 0<br /> sdate = datetime.datetime(year, month, day, hour)<br /> print sdate<br /> #Set directory<br /> dataDir = 'U:/data/cuace_dust/dust_example/2014_case'<br /> outDir = dataDir<br /> infn = os.path.join(dataDir, 'dust_post_'+ sdate.strftime('%Y%m%d%H') + '.ctl')<br /> outfn = os.path.join(dataDir, 'WMO_SDS-WAS_Asian_Center_Model_Forecasting_CUACE-Dust_CMA_' \<br />     + sdate.strftime('%Y-%m-%d') + '.nc')<br /> #Set output X/Y coordinates and projection<br /> toProjInfo = KnownCoordinateSystems.geographic.world.WGS1984<br /> sx = 70.0<br /> xnum = 161<br /> sy = 20<br /> ynum = 71<br /> delta = 0.5<br /> xlist = []<br /> ylist = []<br /> for i in range(0, xnum):<br />     xlist.append(sx + delta * i)<br /> for i in range(0, ynum):<br />     ylist.append(sy + delta * i)<br /> X = jarray.array(xlist, 'd')<br /> Y = jarray.array(ylist, 'd')<br /> <br /> #Read GrADS data file<br /> print 'Open GrADS data file...'<br /> mdi = MeteoDataInfo()<br /> mdi.openGrADSData(infn)<br /> dataInfo = mdi.getDataInfo()<br /> dataInfo.setBigEndian(True)<br /> fromProjInfo = mdi.getProjectionInfo()<br /> tnum = dataInfo.getTimeNum()<br /> mvalue = dataInfo.getMissingValue()<br /> <br /> #Set output nc data file<br /> print 'Create output NC file: ' + outfn<br /> ncfile = NetcdfFileWriter.createNew(NetcdfFileWriter.Version.netcdf3, outfn)<br /> #Add dimensions<br /> print 'Add dimensions...'<br /> xDim = ncfile.addDimension(None, 'lon', xnum)<br /> yDim = ncfile.addDimension(None, 'lat', ynum)<br /> tDim = ncfile.addDimension(None, 'time', tnum)<br /> <br /> #Add global attributes<br /> print 'Add global attributes...'<br /> ncfile.addGroupAttribute(None, Attribute('Conventions', 'CF-1.6'))<br /> ncfile.addGroupAttribute(None, Attribute('Title', 'Sand and dust storm forecasting'))<br /> ncfile.addGroupAttribute(None, Attribute('Model', 'CUACE/Dust'))<br /> ncfile.addGroupAttribute(None, Attribute('Center', 'WMO SDS-WAS Asian Center'))<br /> ncfile.addGroupAttribute(None, Attribute('Agency', 'China Meteorological Administration'))<br /> <br /> #Add variables<br /> xvar = ncfile.addVariable(None, 'lon', DataType.FLOAT, [xDim])<br /> xvar.addAttribute(Attribute('units', 'degrees_east'))<br /> xvar.addAttribute(Attribute('long_name', 'Longitude'))<br /> xvar.addAttribute(Attribute('standard_name', 'longitude'))<br /> xvar.addAttribute(Attribute('axis', 'X'))<br /> yvar = ncfile.addVariable(None, 'lat', DataType.FLOAT, [yDim])<br /> yvar.addAttribute(Attribute('units', 'degrees_north'))<br /> yvar.addAttribute(Attribute('long_name', 'Latitude'))<br /> yvar.addAttribute(Attribute('standard_name', 'latitude'))<br /> yvar.addAttribute(Attribute('axis', 'Y'))<br /> tvar = ncfile.addVariable(None, 'time', DataType.INT, [tDim])<br /> tvar.addAttribute(Attribute('units', 'hours since 1900-01-01 00:00:0.0'))<br /> tvar.addAttribute(Attribute('long_name', 'Time'))<br /> tvar.addAttribute(Attribute('standart_name', 'time'))<br /> tvar.addAttribute(Attribute('axis', 'T'))<br /> #Data variables<br /> vnames = ['load','con','dry','wet','aod']<br /> varlist = []<br /> var = ncfile.addVariable(None, 'LOAD_DUST', DataType.FLOAT, [tDim, yDim, xDim])<br /> var.addAttribute(Attribute('long_name', 'Dust load'))<br /> var.addAttribute(Attribute('units', 'kg/m2'))<br /> var.addAttribute(Attribute('missing_value', -9999.0))<br /> varlist.append(var)<br /> var = ncfile.addVariable(None, 'SCONC_DUST', DataType.FLOAT, [tDim, yDim, xDim])<br /> var.addAttribute(Attribute('long_name', 'Surface dust concentration'))<br /> var.addAttribute(Attribute('units', 'ug/m3'))<br /> var.addAttribute(Attribute('missing_value', -9999.0))<br /> varlist.append(var)<br /> var = ncfile.addVariable(None, 'DDEPO_DUST', DataType.FLOAT, [tDim, yDim, xDim])<br /> var.addAttribute(Attribute('long_name', '3-hour accumulated dry deposition'))<br /> var.addAttribute(Attribute('units', 'kg/m2'))<br /> var.addAttribute(Attribute('missing_value', -9999.0))<br /> varlist.append(var)<br /> var = ncfile.addVariable(None, 'WDEPO_DUST', DataType.FLOAT, [tDim, yDim, xDim])<br /> var.addAttribute(Attribute('long_name', '3-hour accumulated wet deposition'))<br /> var.addAttribute(Attribute('units', 'kg/m2'))<br /> var.addAttribute(Attribute('missing_value', -9999.0))<br /> varlist.append(var)<br /> var = ncfile.addVariable(None, 'AOD550_DUST', DataType.FLOAT, [tDim, yDim, xDim])<br /> var.addAttribute(Attribute('long_name', 'Dust optical depth at 550nm'))<br /> var.addAttribute(Attribute('units', '-'))<br /> var.addAttribute(Attribute('missing_value', -9999.0))<br /> varlist.append(var)<br /> <br /> #Write nc file<br /> ncfile.create()<br /> #Write x,y,z,t variables<br /> print 'Write x variable...'<br /> shape = jarray.array([xnum], 'i')<br /> data = Array.factory(DataType.FLOAT, shape)<br /> for i in range(0,xnum):<br />     data.set(i, X<i>)<br /> ncfile.write(xvar, data)<br /> <br /> print 'Write y variable...'<br /> shape = jarray.array([ynum], 'i')<br /> data = Array.factory(DataType.FLOAT, shape)<br /> for i in range(0,ynum):<br />     data.set(i, Y<i>)<br /> ncfile.write(yvar, data)<br /> <br /> print 'Write time variable...'<br /> format = SimpleDateFormat('yyyy-MM-dd')<br /> sdate = format.parse('1900-01-01')<br /> tvalues = dataInfo.getTimeValues(sdate, 'hours')<br /> shape = jarray.array([tnum], 'i')<br /> data = Array.factory(DataType.INT, shape)<br /> for i in range(0,tnum):<br />     data.set(i, tvalues<i>)<br /> ncfile.write(tvar, data)<br /> <br /> #Write data variables<br /> print 'Write data variable...'<br /> for vname, var in zip(vnames, varlist):<br />     for t in range(0, tnum):<br />         print 'Time: ' + str(t + 1)<br />         mdi.setTimeIndex(t)<br />         gData = mdi.getGridData(vname)<br />         ngData = gData.project(fromProjInfo, toProjInfo, X, Y, ResampleMethods.Bilinear)<br />         origin = jarray.array([t, 0, 0], 'i')<br />         ncfile.write(var, origin, NetCDFDataInfo.gridToArray3D(ngData))<br /> <br /> #Close nc file<br /> ncfile.flush()<br /> ncfile.close()<br /> <br /> print 'Finished'

MeteoInfo

上面转换的netCDF文件绘制模式结果和地面天气现象观测叠加动画图的示例脚本：

1. # coding=utf-8
   
2. #-----------------------------------------------------
   
3. # Author: Yaqiang Wang
   
4. # Date: 2015-3-13
   
5. # Purpose: Read CUACE/Dust netCDF data and MICAPS observation data to plot figures
   
6. # Note: Sample
   
7. #-----------------------------------------------------
   
8. print 'Loading classes...'
   
9. from org.meteoinfo.layout import MapLayout
   
10. from org.meteoinfo.data import GridData
    
11. from org.meteoinfo.data.meteodata import MeteoDataInfo, DrawMeteoData
    
12. from org.meteoinfo.legend import LegendScheme
    
13. from org.meteoinfo.shape import ShapeTypes
    
14. from org.meteoinfo.global.image import AnimatedGifEncoder
    
15. import os.path
    
16. import jarray
    
17. import datetime
    
18. import sys
    
19. from java.util import Date, Calendar, Locale
    
20. from java.text import SimpleDateFormat
    
21. from java.awt import Color
    
22. 
    
23. #Set date
    
24. year = 2013
    
25. month = 3
    
26. day = 1
    
27. hour = 0
    
28. sdate = datetime.datetime(year, month, day, hour)
    
29. 
    
30. #sdate = datetime.date.today()
    
31. #if len(sys.argv) >= 2:
    
32. #    sdate = sdate - datetime.timedelta(days=int(sys.argv[1]))
    
33. #    sdate = sdate + datetime.timedelta(days=1)
    
34. print sdate
    
35. dformat = SimpleDateFormat('HH dd MMM yyy', Locale.ENGLISH)
    
36. dformat1 = SimpleDateFormat('yyMMddHH')
    
37. cal = Calendar.getInstance()
    
38. 
    
39. #Set model
    
40. #model = 'CUACE-DUST_CMA'
    
41. model = 'ADAM2_KMA'
    
42. 
    
43. #Set directory
    
44. dataDir = 'D:/Working/2015/International/SDS_Asian_Region_Center/Model_Verification'
    
45. obsDir = 'U:/data/micaps/2014/plot'
    
46. obsDir = 'E:/MetData/micaps'
    
47. runDir = dataDir
    
48. outDir = os.path.join(dataDir, 'figure')
    
49. if not os.path.exists(outDir):
    
50.     os.mkdir(outDir)
    
51. #Set input/output file names
    
52. infn = os.path.join(dataDir, 'WMO_SDS-WAS_Asian_Center_Model_Forecasting_' + model + '_' \
    
53.     + sdate.strftime('%Y-%m-%d') + '.nc')
    
54. projfn = os.path.join(runDir, 'sds_asian.mip')
    
55. 
    
56. #Plot data
    
57. print 'Plot data...'
    
58. mapLayout = MapLayout()
    
59. mapLayout.loadProjectFile(projfn)
    
60. mf = mapLayout.getActiveMapFrame()
    
61. title = mapLayout.getTexts().get(2)
    
62. legend = mapLayout.getLegends()[0]
    
63. 
    
64. #---- Set weather list - sand and dust storm
    
65. weathers = [6, 7, 8, 9, 30, 31, 32, 33, 34, 35]
    
66. #---- Set weather list - sand and dust storm and haze
    
67. #weathers = [5, 6, 7, 8, 9, 30, 31, 32, 33, 34, 35]
    
68. 
    
69. #---- Create MeteoDataInfo object
    
70. mdi = MeteoDataInfo()
    
71. omdi = MeteoDataInfo()
    
72. 
    
73. #---- Plot loop
    
74. mdi.openNetCDFData(infn)
    
75. lsfn = os.path.join(runDir,'dust_conc.lgs')
    
76. print 'Read data file: ' + infn
    
77. aLS = LegendScheme(ShapeTypes.Polygon)
    
78. aLS.importFromXMLFile(lsfn)
    
79. tnum = mdi.getDataInfo().getTimeNum()
    
80. #tnum = 3
    
81. s = 'SCONC_DUST'
    
82. giffn = os.path.join(outDir, 'V_' + s + '_' + model + '_' + sdate.strftime('%Y%m%d') + '--loop-.gif')
    
83. print giffn
    
84. encoder = AnimatedGifEncoder()
    
85. encoder.setRepeat(0)
    
86. encoder.setDelay(1000)
    
87. encoder.start(giffn)
    
88. sTime = mdi.getDataInfo().getTimes().get(0)
    
89. for t in range(1, tnum):
    
90.     mdi.setTimeIndex(t)
    
91.     aTime = mdi.getDataInfo().getTimes().get(t)
    
92.     cal.setTime(aTime)
    
93.     cal.add(Calendar.HOUR, 8)
    
94.     bjTime = cal.getTime()
    
95.     #---- Open observation weather data   
    
96.     obsfn = os.path.join(obsDir, dformat1.format(bjTime) + '.000')
    
97.     print obsfn
    
98.     if not os.path.exists(obsfn):
    
99.         continue
    
100.     omdi.openMICAPSData(obsfn)
     
101.     wData = omdi.getStationData('WeatherNow')
     
102.     weatherLayer = DrawMeteoData.createWeatherSymbolLayer(wData, weathers, 'Weather')
     
103.     #for lb in weatherLayer.getLegendScheme().getLegendBreaks():
     
104.     #    lb.setColor(Color.red)
     
105.     weatherLayer.setAvoidCollision(False)
     
106.     mf.removeMeteoLayers()
     
107.     mf.addLayer(weatherLayer)
     
108.     #---- Get grid data and create a shaded layer
     
109.     gData = mdi.getGridData(s)  
     
110.     aLayer = DrawMeteoData.createShadedLayer(gData, aLS, 'Forecasting_' + s, 'Data', True)
     
111.     aLayer.setProjInfo(mdi.getProjectionInfo())   
     
112.     mf.addLayer(aLayer)
     
113.     mf.moveLayer(aLayer, 0)
     
114.     #---- Set title
     
115.     title.setLabelText('Run: ' + dformat.format(sTime) + '    Valid: ' + dformat.format(aTime) \
     
116.         + '(H+' + str(t * 3) + ')')
     
117.     #---- Set legend
     
118.     legend.setLegendLayer(aLayer)
     
119.     mapLayout.paintGraphics()
     
120.     encoder.addFrame(mapLayout.getViewImage())
     
121.     figurefn = os.path.join(outDir, 'V_' + model + '_' + s + '_' + dformat1.format(aTime) + '.png')
     
122.     print 'Output figure: ' + figurefn
     
123.     mapLayout.exportToPicture(figurefn)
     
124. 
     
125. encoder.finish()
     
126. print 'Finished'

# coding=utf-8<br /> #-----------------------------------------------------<br /> # Author: Yaqiang Wang<br /> # Date: 2015-3-13<br /> # Purpose: Read CUACE/Dust netCDF data and MICAPS observation data to plot figures<br /> # Note: Sample<br /> #-----------------------------------------------------<br /> print 'Loading classes...'<br /> from org.meteoinfo.layout import MapLayout<br /> from org.meteoinfo.data import GridData<br /> from org.meteoinfo.data.meteodata import MeteoDataInfo, DrawMeteoData<br /> from org.meteoinfo.legend import LegendScheme<br /> from org.meteoinfo.shape import ShapeTypes<br /> from org.meteoinfo.global.image import AnimatedGifEncoder<br /> import os.path<br /> import jarray<br /> import datetime<br /> import sys<br /> from java.util import Date, Calendar, Locale<br /> from java.text import SimpleDateFormat<br /> from java.awt import Color<br /> <br /> #Set date<br /> year = 2013<br /> month = 3<br /> day = 1<br /> hour = 0<br /> sdate = datetime.datetime(year, month, day, hour)<br /> <br /> #sdate = datetime.date.today()<br /> #if len(sys.argv) >= 2:<br /> #    sdate = sdate - datetime.timedelta(days=int(sys.argv[1]))<br /> #    sdate = sdate + datetime.timedelta(days=1)<br /> print sdate<br /> dformat = SimpleDateFormat('HH dd MMM yyy', Locale.ENGLISH)<br /> dformat1 = SimpleDateFormat('yyMMddHH')<br /> cal = Calendar.getInstance()<br /> <br /> #Set model<br /> #model = 'CUACE-DUST_CMA'<br /> model = 'ADAM2_KMA'<br /> <br /> #Set directory<br /> dataDir = 'D:/Working/2015/International/SDS_Asian_Region_Center/Model_Verification'<br /> obsDir = 'U:/data/micaps/2014/plot'<br /> obsDir = 'E:/MetData/micaps'<br /> runDir = dataDir<br /> outDir = os.path.join(dataDir, 'figure')<br /> if not os.path.exists(outDir):<br />     os.mkdir(outDir)<br /> #Set input/output file names<br /> infn = os.path.join(dataDir, 'WMO_SDS-WAS_Asian_Center_Model_Forecasting_' + model + '_' \<br />     + sdate.strftime('%Y-%m-%d') + '.nc')<br /> projfn = os.path.join(runDir, 'sds_asian.mip')<br /> <br /> #Plot data<br /> print 'Plot data...'<br /> mapLayout = MapLayout()<br /> mapLayout.loadProjectFile(projfn)<br /> mf = mapLayout.getActiveMapFrame()<br /> title = mapLayout.getTexts().get(2)<br /> legend = mapLayout.getLegends()[0]<br /> <br /> #---- Set weather list - sand and dust storm<br /> weathers = [6, 7, 8, 9, 30, 31, 32, 33, 34, 35]<br /> #---- Set weather list - sand and dust storm and haze<br /> #weathers = [5, 6, 7, 8, 9, 30, 31, 32, 33, 34, 35]<br /> <br /> #---- Create MeteoDataInfo object<br /> mdi = MeteoDataInfo()<br /> omdi = MeteoDataInfo()<br /> <br /> #---- Plot loop<br /> mdi.openNetCDFData(infn)<br /> lsfn = os.path.join(runDir,'dust_conc.lgs')<br /> print 'Read data file: ' + infn<br /> aLS = LegendScheme(ShapeTypes.Polygon)<br /> aLS.importFromXMLFile(lsfn)<br /> tnum = mdi.getDataInfo().getTimeNum()<br /> #tnum = 3<br /> s = 'SCONC_DUST'<br /> giffn = os.path.join(outDir, 'V_' + s + '_' + model + '_' + sdate.strftime('%Y%m%d') + '--loop-.gif')<br /> print giffn<br /> encoder = AnimatedGifEncoder()<br /> encoder.setRepeat(0)<br /> encoder.setDelay(1000)<br /> encoder.start(giffn)<br /> sTime = mdi.getDataInfo().getTimes().get(0)<br /> for t in range(1, tnum):<br />     mdi.setTimeIndex(t)<br />     aTime = mdi.getDataInfo().getTimes().get(t)<br />     cal.setTime(aTime)<br />     cal.add(Calendar.HOUR, 8)<br />     bjTime = cal.getTime()<br />     #---- Open observation weather data    <br />     obsfn = os.path.join(obsDir, dformat1.format(bjTime) + '.000')<br />     print obsfn<br />     if not os.path.exists(obsfn):<br />         continue<br />     omdi.openMICAPSData(obsfn)<br />     wData = omdi.getStationData('WeatherNow')<br />     weatherLayer = DrawMeteoData.createWeatherSymbolLayer(wData, weathers, 'Weather')<br />     #for lb in weatherLayer.getLegendScheme().getLegendBreaks():<br />     #    lb.setColor(Color.red)<br />     weatherLayer.setAvoidCollision(False)<br />     mf.removeMeteoLayers()<br />     mf.addLayer(weatherLayer)<br />     #---- Get grid data and create a shaded layer<br />     gData = mdi.getGridData(s)  <br />     aLayer = DrawMeteoData.createShadedLayer(gData, aLS, 'Forecasting_' + s, 'Data', True)<br />     aLayer.setProjInfo(mdi.getProjectionInfo())    <br />     mf.addLayer(aLayer)<br />     mf.moveLayer(aLayer, 0)<br />     #---- Set title<br />     title.setLabelText('Run: ' + dformat.format(sTime) + '    Valid: ' + dformat.format(aTime) \<br />         + '(H+' + str(t * 3) + ')')<br />     #---- Set legend<br />     legend.setLegendLayer(aLayer)<br />     mapLayout.paintGraphics()<br />     encoder.addFrame(mapLayout.getViewImage())<br />     figurefn = os.path.join(outDir, 'V_' + model + '_' + s + '_' + dformat1.format(aTime) + '.png')<br />     print 'Output figure: ' + figurefn<br />     mapLayout.exportToPicture(figurefn)<br /> <br /> encoder.finish()<br /> print 'Finished'

meteomm

谢谢LZ，我研究下。

mzy1027

这个帖子怎么收藏？？？

lwswl

知道了。。。

lwswl

mzy1027

怎么修改用户名？

精卫

大大的点个赞

人生路

MeteoInfo 发表于 2015-4-16 17:39
上面转换的netCDF文件绘制模式结果和地面天气现象观测叠加动画图的示例脚本：

王老师，这个meteoinfo能用地图文件裁剪nc数据得到底图的区域吗

MeteoInfo

人生路 发表于 2015-5-3 12:25
王老师，这个meteoinfo能用地图文件裁剪nc数据得到底图的区域吗

可以