MeteoInfo脚本示例：GrADS to netCDF

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

#-----------------------------------------------------
# Author: Yaqiang Wang
# Date: 2015-3-12
# Purpose: Convert CUACE grads data to netCDF (CUACE/Dust)
# Note: Sample
#-----------------------------------------------------
print 'Loading classes...'
from org.meteoinfo.data import GridData
from org.meteoinfo.data import DataMath
from org.meteoinfo.data.meteodata import MeteoDataInfo
from org.meteoinfo.geoprocess.analysis import ResampleMethods
from org.meteoinfo.data.meteodata.netcdf import NetCDFDataInfo
from org.meteoinfo.projection import ProjectionInfo
from org.meteoinfo.projection import ProjectionNames
from org.meteoinfo.projection import KnownCoordinateSystems
from ucar.nc2 import NetcdfFileWriter
from ucar.nc2 import Attribute
from ucar.ma2 import DataType
from ucar.ma2 import Array
import os.path
import jarray
import datetime
from java.util import Date
from java.text import SimpleDateFormat

#Set date
year = 2014
month = 4
day = 23
hour = 0
sdate = datetime.datetime(year, month, day, hour)
print sdate
#Set directory
dataDir = 'U:/data/cuace_dust/dust_example/2014_case'
outDir = dataDir
infn = os.path.join(dataDir, 'dust_post_'+ sdate.strftime('%Y%m%d%H') + '.ctl')
outfn = os.path.join(dataDir, 'WMO_SDS-WAS_Asian_Center_Model_Forecasting_CUACE-Dust_CMA_' 
    + sdate.strftime('%Y-%m-%d') + '.nc')
#Set output X/Y coordinates and projection
toProjInfo = KnownCoordinateSystems.geographic.world.WGS1984
sx = 70.0
xnum = 161
sy = 20
ynum = 71
delta = 0.5
xlist = []
ylist = []
for i in range(0, xnum):
    xlist.append(sx + delta * i)
for i in range(0, ynum):
    ylist.append(sy + delta * i)
X = jarray.array(xlist, 'd')
Y = jarray.array(ylist, 'd')

#Read GrADS data file
print 'Open GrADS data file...'
mdi = MeteoDataInfo()
mdi.openGrADSData(infn)
dataInfo = mdi.getDataInfo()
dataInfo.setBigEndian(True)
fromProjInfo = mdi.getProjectionInfo()
tnum = dataInfo.getTimeNum()
mvalue = dataInfo.getMissingValue()

#Set output nc data file
print 'Create output NC file: ' + outfn
ncfile = NetcdfFileWriter.createNew(NetcdfFileWriter.Version.netcdf3, outfn)
#Add dimensions
print 'Add dimensions...'
xDim = ncfile.addDimension(None, 'lon', xnum)
yDim = ncfile.addDimension(None, 'lat', ynum)
tDim = ncfile.addDimension(None, 'time', tnum)

#Add global attributes
print 'Add global attributes...'
ncfile.addGroupAttribute(None, Attribute('Conventions', 'CF-1.6'))
ncfile.addGroupAttribute(None, Attribute('Title', 'Sand and dust storm forecasting'))
ncfile.addGroupAttribute(None, Attribute('Model', 'CUACE/Dust'))
ncfile.addGroupAttribute(None, Attribute('Center', 'WMO SDS-WAS Asian Center'))
ncfile.addGroupAttribute(None, Attribute('Agency', 'China Meteorological Administration'))

#Add variables
xvar = ncfile.addVariable(None, 'lon', DataType.FLOAT, [xDim])
xvar.addAttribute(Attribute('units', 'degrees_east'))
xvar.addAttribute(Attribute('long_name', 'Longitude'))
xvar.addAttribute(Attribute('standard_name', 'longitude'))
xvar.addAttribute(Attribute('axis', 'X'))
yvar = ncfile.addVariable(None, 'lat', DataType.FLOAT, [yDim])
yvar.addAttribute(Attribute('units', 'degrees_north'))
yvar.addAttribute(Attribute('long_name', 'Latitude'))
yvar.addAttribute(Attribute('standard_name', 'latitude'))
yvar.addAttribute(Attribute('axis', 'Y'))
tvar = ncfile.addVariable(None, 'time', DataType.INT, [tDim])
tvar.addAttribute(Attribute('units', 'hours since 1900-01-01 00:00:0.0'))
tvar.addAttribute(Attribute('long_name', 'Time'))
tvar.addAttribute(Attribute('standart_name', 'time'))
tvar.addAttribute(Attribute('axis', 'T'))
#Data variables
vnames = ['load','con','dry','wet','aod']
varlist = []
var = ncfile.addVariable(None, 'LOAD_DUST', DataType.FLOAT, [tDim, yDim, xDim])
var.addAttribute(Attribute('long_name', 'Dust load'))
var.addAttribute(Attribute('units', 'kg/m2'))
var.addAttribute(Attribute('missing_value', -9999.0))
varlist.append(var)
var = ncfile.addVariable(None, 'SCONC_DUST', DataType.FLOAT, [tDim, yDim, xDim])
var.addAttribute(Attribute('long_name', 'Surface dust concentration'))
var.addAttribute(Attribute('units', 'ug/m3'))
var.addAttribute(Attribute('missing_value', -9999.0))
varlist.append(var)
var = ncfile.addVariable(None, 'DDEPO_DUST', DataType.FLOAT, [tDim, yDim, xDim])
var.addAttribute(Attribute('long_name', '3-hour accumulated dry deposition'))
var.addAttribute(Attribute('units', 'kg/m2'))
var.addAttribute(Attribute('missing_value', -9999.0))
varlist.append(var)
var = ncfile.addVariable(None, 'WDEPO_DUST', DataType.FLOAT, [tDim, yDim, xDim])
var.addAttribute(Attribute('long_name', '3-hour accumulated wet deposition'))
var.addAttribute(Attribute('units', 'kg/m2'))
var.addAttribute(Attribute('missing_value', -9999.0))
varlist.append(var)
var = ncfile.addVariable(None, 'AOD550_DUST', DataType.FLOAT, [tDim, yDim, xDim])
var.addAttribute(Attribute('long_name', 'Dust optical depth at 550nm'))
var.addAttribute(Attribute('units', '-'))
var.addAttribute(Attribute('missing_value', -9999.0))
varlist.append(var)

#Write nc file
ncfile.create()
#Write x,y,z,t variables
print 'Write x variable...'
shape = jarray.array([xnum], 'i')
data = Array.factory(DataType.FLOAT, shape)
for i in range(0,xnum):
    data.set(i, X[i])
ncfile.write(xvar, data)

print 'Write y variable...'
shape = jarray.array([ynum], 'i')
data = Array.factory(DataType.FLOAT, shape)
for i in range(0,ynum):
    data.set(i, Y[i])
ncfile.write(yvar, data)

print 'Write time variable...'
format = SimpleDateFormat('yyyy-MM-dd')
sdate = format.parse('1900-01-01')
tvalues = dataInfo.getTimeValues(sdate, 'hours')
shape = jarray.array([tnum], 'i')
data = Array.factory(DataType.INT, shape)
for i in range(0,tnum):
    data.set(i, tvalues[i])
ncfile.write(tvar, data)

#Write data variables
print 'Write data variable...'
for vname, var in zip(vnames, varlist):
    for t in range(0, tnum):
        print 'Time: ' + str(t + 1)
        mdi.setTimeIndex(t)
        gData = mdi.getGridData(vname)
        ngData = gData.project(fromProjInfo, toProjInfo, X, Y, ResampleMethods.Bilinear)
        origin = jarray.array([t, 0, 0], 'i')
        ncfile.write(var, origin, NetCDFDataInfo.gridToArray3D(ngData))

#Close nc file
ncfile.flush()
ncfile.close()

print 'Finished'

　　

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

# coding=utf-8
#-----------------------------------------------------
# Author: Yaqiang Wang
# Date: 2015-3-13
# Purpose: Read CUACE/Dust netCDF data and MICAPS observation data to plot figures
# Note: Sample
#-----------------------------------------------------
print 'Loading classes...'
from org.meteoinfo.layout import MapLayout
from org.meteoinfo.data import GridData
from org.meteoinfo.data.meteodata import MeteoDataInfo, DrawMeteoData
from org.meteoinfo.legend import LegendScheme
from org.meteoinfo.shape import ShapeTypes
from org.meteoinfo.global.image import AnimatedGifEncoder
import os.path
import jarray
import datetime
import sys
from java.util import Date, Calendar, Locale
from java.text import SimpleDateFormat
from java.awt import Color

#Set date
year = 2013
month = 3
day = 1
hour = 0
sdate = datetime.datetime(year, month, day, hour)

#sdate = datetime.date.today()
#if len(sys.argv) >= 2:
#    sdate = sdate - datetime.timedelta(days=int(sys.argv[1]))
#    sdate = sdate + datetime.timedelta(days=1)
print sdate
dformat = SimpleDateFormat('HH dd MMM yyy', Locale.ENGLISH)
dformat1 = SimpleDateFormat('yyMMddHH')
cal = Calendar.getInstance()

#Set model
#model = 'CUACE-DUST_CMA'
model = 'ADAM2_KMA'

#Set directory
dataDir = 'D:/Working/2015/International/SDS_Asian_Region_Center/Model_Verification'
obsDir = 'U:/data/micaps/2014/plot'
obsDir = 'E:/MetData/micaps'
runDir = dataDir
outDir = os.path.join(dataDir, 'figure')
if not os.path.exists(outDir):
    os.mkdir(outDir)
#Set input/output file names
infn = os.path.join(dataDir, 'WMO_SDS-WAS_Asian_Center_Model_Forecasting_' + model + '_' 
    + sdate.strftime('%Y-%m-%d') + '.nc')
projfn = os.path.join(runDir, 'sds_asian.mip')

#Plot data
print 'Plot data...'
mapLayout = MapLayout()
mapLayout.loadProjectFile(projfn)
mf = mapLayout.getActiveMapFrame()
title = mapLayout.getTexts().get(2)
legend = mapLayout.getLegends()[0]

#---- Set weather list - sand and dust storm
weathers = [6, 7, 8, 9, 30, 31, 32, 33, 34, 35]
#---- Set weather list - sand and dust storm and haze
#weathers = [5, 6, 7, 8, 9, 30, 31, 32, 33, 34, 35]

#---- Create MeteoDataInfo object
mdi = MeteoDataInfo()
omdi = MeteoDataInfo()

#---- Plot loop
mdi.openNetCDFData(infn)
lsfn = os.path.join(runDir,'dust_conc.lgs')
print 'Read data file: ' + infn
aLS = LegendScheme(ShapeTypes.Polygon)
aLS.importFromXMLFile(lsfn)
tnum = mdi.getDataInfo().getTimeNum()
#tnum = 3
s = 'SCONC_DUST'
giffn = os.path.join(outDir, 'V_' + s + '_' + model + '_' + sdate.strftime('%Y%m%d') + '--loop-.gif')
print giffn
encoder = AnimatedGifEncoder()
encoder.setRepeat(0)
encoder.setDelay(1000)
encoder.start(giffn)
sTime = mdi.getDataInfo().getTimes().get(0)
for t in range(1, tnum):
    mdi.setTimeIndex(t)
    aTime = mdi.getDataInfo().getTimes().get(t)
    cal.setTime(aTime)
    cal.add(Calendar.HOUR, 8)
    bjTime = cal.getTime()
    #---- Open observation weather data    
    obsfn = os.path.join(obsDir, dformat1.format(bjTime) + '.000')
    print obsfn
    if not os.path.exists(obsfn):
        continue
    omdi.openMICAPSData(obsfn)
    wData = omdi.getStationData('WeatherNow')
    weatherLayer = DrawMeteoData.createWeatherSymbolLayer(wData, weathers, 'Weather')
    #for lb in weatherLayer.getLegendScheme().getLegendBreaks():
    #    lb.setColor(Color.red)
    weatherLayer.setAvoidCollision(False)
    mf.removeMeteoLayers()
    mf.addLayer(weatherLayer)
    #---- Get grid data and create a shaded layer
    gData = mdi.getGridData(s)  
    aLayer = DrawMeteoData.createShadedLayer(gData, aLS, 'Forecasting_' + s, 'Data', True)
    aLayer.setProjInfo(mdi.getProjectionInfo())    
    mf.addLayer(aLayer)
    mf.moveLayer(aLayer, 0)
    #---- Set title
    title.setLabelText('Run: ' + dformat.format(sTime) + '    Valid: ' + dformat.format(aTime) 
        + '(H+' + str(t * 3) + ')')
    #---- Set legend
    legend.setLegendLayer(aLayer)
    mapLayout.paintGraphics()
    encoder.addFrame(mapLayout.getViewImage())
    figurefn = os.path.join(outDir, 'V_' + model + '_' + s + '_' + dformat1.format(aTime) + '.png')
    print 'Output figure: ' + figurefn
    mapLayout.exportToPicture(figurefn)

encoder.finish()
print 'Finished'