MeteoInfo脚本教程（十六）：格点数据运算

MeteoInfo · 发表于 2014-12-21 00:20:21

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

数据分析计算也是MeteoInfo重点考虑的功能，主要体现在格点和站点数据的各种数学运算。这里先讲讲格点数据的运算。

我一直认为Java语言有一个重大缺陷，那就是不支持运算符重载，这造成无法用简单的运算符号（比如+, -, *, /）来写公式进行自定义数据类的运算。所以在GridData类中定义了一些方法来进行四则运算：add, sub, mul, div。还需要知道一个重要的类DataMath（在org.meteoinfo.data包中），利用该类可以进行更复杂的运算，比如指数、幂、三角函数运算等等。这样处理也可以完成运算符重载能做的工作，但是公式写起来就没那么简洁明了了。比如如果有运算符重载功能，可以写下面的公式（air是GridData类的对象）：
es = 6.112*DataMath.exp(17.67*(air-273.16)/(air-29.65))

没有运算符重载功能只能处理成这样：
es = DataMath.exp(air.sub(273.16).mul(17.67).div(air.sub(29.65))).mul(6.112)

这是题外话，Java的开发人员不认为这是一个缺陷，而且也没有加入运算符重载功能的打算，所以公式写得难看些也只能这样了。

下面是一个利用NCEP数据计算水汽通量散度的脚本，包含了比较多的格点运算。

#-----------------------------------------------------
# Author: Yaqiang Wang
# Date: 2014-12-20
# Purpose: Calculate water vapor flux
# Note: Sample
#-----------------------------------------------------
#---- Import classes
print 'Import classes...'
from org.meteoinfo.layout import MapLayout
from org.meteoinfo.data import DataMath
from org.meteoinfo.data.mapdata import MapDataManage
from org.meteoinfo.data.meteodata import MeteoDataInfo, DrawMeteoData
from org.meteoinfo.legend import LegendManage, LegendType
from org.meteoinfo.layout import LegendStyles
from org.meteoinfo.global import Extent
import os.path
from java.text import SimpleDateFormat
from java.awt import Color
from javax.swing import JFrame
#---- Set directories
print 'Set directories...'
baseDir = 'D:/MyProgram/Distribution/java/MeteoInfo/MeteoInfo'
dataDir = 'D:/Temp/nc'
mapDir = os.path.join(baseDir, 'map')
#---- Create MapLayout object
mapLayout = MapLayout()
mapFrame = mapLayout.getActiveMapFrame()
layoutMap = mapLayout.getActiveLayoutMap()
#---- Load country layer
print 'Load country layer...'
countryLayer = MapDataManage.loadLayer(os.path.join(mapDir, 'country1.shp'))
lb = countryLayer.getLegendScheme().getLegendBreaks().get(0)
lb.setDrawFill(False)
lb.setOutlineColor(Color.blue)
mapFrame.addLayer(countryLayer)
#---- Open netCDF data files
print 'Open netCDF data files...'
dataAir = MeteoDataInfo()
dataUwnd = MeteoDataInfo()
dataVwnd = MeteoDataInfo()
dataRhum = MeteoDataInfo()
dataAir.openNetCDFData(os.path.join(dataDir, 'air.2011.nc'))
dataUwnd.openNetCDFData(os.path.join(dataDir, 'uwnd.2011.nc'))
dataVwnd.openNetCDFData(os.path.join(dataDir, 'vwnd.2011.nc'))
dataRhum.openNetCDFData(os.path.join(dataDir, 'rhum.2011.nc'))
#---- Calculation
#---- Set date to June 23
#tIdx = 171
tIdx = 173
dataAir.setTimeIndex(tIdx);
dataUwnd.setTimeIndex(tIdx);
dataVwnd.setTimeIndex(tIdx);
dataRhum.setTimeIndex(tIdx);
#---- Set level to 700hPa
lIdx = 3
dataAir.setLevelIndex(lIdx);
dataUwnd.setLevelIndex(lIdx);
dataVwnd.setLevelIndex(lIdx);
dataRhum.setLevelIndex(lIdx);
#---- Get grid data
print 'Get grid data...'
air = dataAir.getGridData('air')
uwnd = dataUwnd.getGridData('uwnd')
vwnd = dataVwnd.getGridData('vwnd')
rhum = dataRhum.getGridData('rhum')
#---- Calculate
print 'Calculation...'
prs = 700
g = 9.8
es = DataMath.exp(air.sub(273.16).mul(17.67).div(air.sub(29.65))).mul(6.112)
#es = 6.112*DataMath.exp(17.67*(air-273.16)/(air-29.65))
qs = es.mul(0.62197).div(DataMath.sub(prs, es.mul(0.378)))
#qs = 0.62197*es/(prs-0.378*es)
q = qs.mul(rhum).div(100)
#q = qs*rhum/100
qhdivg = DataMath.hdivg(q.mul(uwnd).div(g), q.mul(vwnd).div(g))
#qhdivg = DataMath.Hdivg(q*uwnd/g,q*vwnd/g)
qv = rhum.mul(es).div(100)
#qv = rhum*es/100
uv = DataMath.magnitude(uwnd, vwnd)
#uv = DataMath.Magnitude(uwnd, vwnd)
uvq = uv.mul(qv).div(9.8*1000)
#uvq = uv*qv/(9.8*1000)
#---- Create data layer
print 'Create data layer...'
dataLayer = DrawMeteoData.createShadedLayer(qhdivg, "WaterVaporFlux", "Flux", False)
#---- Add layer
print 'Add layers...'
mapFrame.addLayer(dataLayer)
mapFrame.moveLayer(dataLayer, 0)
#---- Zoom data
mapLayout.getActiveLayoutMap().zoomToExtentLonLatEx(Extent(0,360,-90.1,90.1))
#---- Set MapLayout
format = SimpleDateFormat('yyyy-MM-dd')
aTime = dataAir.getDataInfo().getTimes().get(tIdx)
mapLayout.addText('Water Vapor Flux Divergence - ' + format.format(aTime), 320, 30, 'Arial', 16)
aLegend = mapLayout.addLegend(650, 100)
aLegend.setLegendStyle(LegendStyles.Bar_Vertical)
aLegend.setLegendLayer(dataLayer)
layoutMap.setGridXDelt(60)
layoutMap.setGridYDelt(30)
layoutMap.setDrawGridLine(False)
mapLayout.paintGraphics()
frame = JFrame('MeteoInfo Script Sample', size = (750, 530))
frame.add(mapLayout)
frame.visible = True
print 'Finished!'

#----------------------------------------------------- 
# Author: Yaqiang Wang 
# Date: 2014-12-20 
# Purpose: Calculate water vapor flux 
# Note: Sample 
#----------------------------------------------------- 
#---- Import classes 
print 'Import classes...' 
from org.meteoinfo.layout import MapLayout 
from org.meteoinfo.data import DataMath 
from org.meteoinfo.data.mapdata import MapDataManage 
from org.meteoinfo.data.meteodata import MeteoDataInfo, DrawMeteoData 
from org.meteoinfo.legend import LegendManage, LegendType 
from org.meteoinfo.layout import LegendStyles 
from org.meteoinfo.global import Extent 
import os.path 
from java.text import SimpleDateFormat 
from java.awt import Color 
from javax.swing import JFrame 
 
#---- Set directories 
print 'Set directories...' 
baseDir = 'D:/MyProgram/Distribution/java/MeteoInfo/MeteoInfo' 
dataDir = 'D:/Temp/nc' 
mapDir = os.path.join(baseDir, 'map') 
 
#---- Create MapLayout object 
mapLayout = MapLayout() 
mapFrame = mapLayout.getActiveMapFrame() 
layoutMap = mapLayout.getActiveLayoutMap() 
 
#---- Load country layer 
print 'Load country layer...' 
countryLayer = MapDataManage.loadLayer(os.path.join(mapDir, 'country1.shp')) 
lb = countryLayer.getLegendScheme().getLegendBreaks().get(0) 
lb.setDrawFill(False) 
lb.setOutlineColor(Color.blue) 
mapFrame.addLayer(countryLayer) 
 
#---- Open netCDF data files 
print 'Open netCDF data files...' 
dataAir = MeteoDataInfo() 
dataUwnd = MeteoDataInfo() 
dataVwnd = MeteoDataInfo() 
dataRhum = MeteoDataInfo() 
dataAir.openNetCDFData(os.path.join(dataDir, 'air.2011.nc')) 
dataUwnd.openNetCDFData(os.path.join(dataDir, 'uwnd.2011.nc')) 
dataVwnd.openNetCDFData(os.path.join(dataDir, 'vwnd.2011.nc')) 
dataRhum.openNetCDFData(os.path.join(dataDir, 'rhum.2011.nc')) 
 
#---- Calculation 
#---- Set date to June 23 
#tIdx = 171 
tIdx = 173 
dataAir.setTimeIndex(tIdx); 
dataUwnd.setTimeIndex(tIdx); 
dataVwnd.setTimeIndex(tIdx); 
dataRhum.setTimeIndex(tIdx); 
 
#---- Set level to 700hPa 
lIdx = 3 
dataAir.setLevelIndex(lIdx); 
dataUwnd.setLevelIndex(lIdx); 
dataVwnd.setLevelIndex(lIdx); 
dataRhum.setLevelIndex(lIdx); 
 
#---- Get grid data 
print 'Get grid data...' 
air = dataAir.getGridData('air') 
uwnd = dataUwnd.getGridData('uwnd') 
vwnd = dataVwnd.getGridData('vwnd') 
rhum = dataRhum.getGridData('rhum') 
 
#---- Calculate 
print 'Calculation...' 
prs = 700 
g = 9.8 
es = DataMath.exp(air.sub(273.16).mul(17.67).div(air.sub(29.65))).mul(6.112) 
#es = 6.112*DataMath.exp(17.67*(air-273.16)/(air-29.65)) 
qs = es.mul(0.62197).div(DataMath.sub(prs, es.mul(0.378))) 
#qs = 0.62197*es/(prs-0.378*es) 
q = qs.mul(rhum).div(100) 
#q = qs*rhum/100 
qhdivg = DataMath.hdivg(q.mul(uwnd).div(g), q.mul(vwnd).div(g)) 
#qhdivg = DataMath.Hdivg(q*uwnd/g,q*vwnd/g) 
qv = rhum.mul(es).div(100) 
#qv = rhum*es/100 
uv = DataMath.magnitude(uwnd, vwnd) 
#uv = DataMath.Magnitude(uwnd, vwnd) 
uvq = uv.mul(qv).div(9.8*1000) 
#uvq = uv*qv/(9.8*1000) 
 
#---- Create data layer 
print 'Create data layer...' 
dataLayer = DrawMeteoData.createShadedLayer(qhdivg, "WaterVaporFlux", "Flux", False) 
 
#---- Add layer 
print 'Add layers...' 
mapFrame.addLayer(dataLayer) 
mapFrame.moveLayer(dataLayer, 0) 
#---- Zoom data 
mapLayout.getActiveLayoutMap().zoomToExtentLonLatEx(Extent(0,360,-90.1,90.1)) 
#---- Set MapLayout 
format = SimpleDateFormat('yyyy-MM-dd') 
aTime = dataAir.getDataInfo().getTimes().get(tIdx) 
mapLayout.addText('Water Vapor Flux Divergence - ' + format.format(aTime), 320, 30, 'Arial', 16) 
aLegend = mapLayout.addLegend(650, 100) 
aLegend.setLegendStyle(LegendStyles.Bar_Vertical) 
aLegend.setLegendLayer(dataLayer) 
layoutMap.setGridXDelt(60) 
layoutMap.setGridYDelt(30) 
layoutMap.setDrawGridLine(False) 
mapLayout.paintGraphics() 
 
frame = JFrame('MeteoInfo Script Sample', size = (750, 530)) 
frame.add(mapLayout) 
frame.visible = True 
print 'Finished!'

远山 · 发表于 2014-12-21 09:18:34

呵呵，王老师

兰溪之水 · 发表于 2014-12-21 21:54:05

王老师，您的软件是否能后台事先设置好一些默认设置，然后前台写的代码能尽量的简洁？
如果能做成像GrADS那样能用简单的脚本就可以绘制常用的图，我想王老师的软件肯定会更受欢迎。

MeteoInfo · 发表于 2014-12-21 22:02:22

兰溪之水发表于 2014-12-21 21:54
王老师，您的软件是否能后台事先设置好一些默认设置，然后前台写的代码能尽量的简洁？
如果能做成像GrADS ...

最早C#版MeteoInfo专门设计了一个MIApp类就是想封装一些常用功能，让脚本更简单些。不过这也是双刃剑，首先牺牲了灵活性，再者不利于用户真正了解MeteoInfo库的结构和功能。

当然，可能很多人不愿意学习复杂的东西，以后也会考虑在Java版中做这样的常用功能封装类。

		自动登录	找回密码
密码			立即注册

MeteoInfo脚本教程（十六）：格点数据运算

登录后查看更多精彩内容~

本帖被以下淘专辑推荐:

相关帖子