MeteoInfo脚本示例:GrADS to netCDF
阅读原文时间:2023年07月10日阅读:4

这里给出一个将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'