MeteoInfo脚本示例 - 从再分析数据中提取站点数据

MeteoInfo · 发表于 2014-9-14 19:21:38

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

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

x

本帖最后由 MeteoInfo 于 2014-9-14 19:23 编辑

从格点数据中提取某个点的数据其实就是一个插值过程，通常用双线性插值，MeteoInfo的格点数据类GridData有toStation(double x, double y)方法用双线性插值方法从格点数据中获取某个点的数据。下面两个脚本程序演示了从GDAS再分析数据中获取某点的边界层高度和三维风向风速数据，示例程序用的数据格式是ARL（HYSPLIT模式用的），不过只要是MeteoInfo支持的数据格式（比如netCDF, GRIB等）都可以用下面的脚本程序（只需做简单的改动）。脚本程序里有比较丰富的注释，这里就不多说了。需要说明在Java的Calendar中月份是从0开始的，因此月设为10其实是11月。

#--------------------------------------------------------
# Author: Yaqiang Wang
# Date: 2014-9-11
# Purpose: Extract station PBL data from ARL meteorological data
# Note: Sample
#-----------------------------------------------------------
from org.meteoinfo.data.meteodata import MeteoDataInfo
import os.path
from java.util import Date
from java.util import Calendar
from java.util import Locale
from java.text import SimpleDateFormat
#---- Set directories
dataDir = 'P:/MetData/MetData/2011/'
#---- Set times - this case show seasonal variation
cal = Calendar.getInstance()
cal.set(2011, 10, 1)
sTime = cal.getTime()
cal.set(2011, 11, 1)
eTime = cal.getTime()
aTime = sTime
format = SimpleDateFormat('MMMyy', Locale.ENGLISH)
format1 = SimpleDateFormat('yyyy-MM-dd HH:00')
format2 = SimpleDateFormat('yyyyMMddHH')
cal.setTime(sTime)
#---- Create a MeteoDataInfo object
mydata = MeteoDataInfo()
#---- Set lon/lat of the location
lon = 20.5
lat = 40.5
#---- Set output text file
outpath = os.path.join(dataDir, 'PBL_data-1.txt')
outf = open(outpath, 'w')
outf.write('Time,PBLH')
#---- Loop
i = 0
sum = 0.0
while aTime <= eTime:
for w in range(5):
inFile = dataDir + "gdas1." + format.format(aTime) + ".w" + str(w + 1)
print inFile
if os.path.isfile(inFile):
mydata.openARLData(inFile) #Open ARL data file
tNum = mydata.getDataInfo().getTimeNum()
print tNum
for t in range(tNum):
dTime = mydata.getDataInfo().getTimes().get(t)
print format1.format(dTime)
mydata.setTimeIndex(t) #Set time index
gData = mydata.getGridData("PBLH") #Get PBLH grid data
pblh = gData.toStation(lon, lat) #Intepolate PBLH data to the location
print 'PBLH = %10.2f' %pblh
outf.write(os.linesep)
outf.write(format2.format(dTime) + ',%-7.2f ' %pblh)
sum += pblh
i += 1
cal.add(Calendar.MONTH, 1)
aTime = cal.getTime()
outf.close()
#---- Calculate seasonal average PBLH value
ave = sum / i
print 'Average PBLH: %10.2f' %ave
print 'Finished!'

#--------------------------------------------------------         
# Author: Yaqiang Wang                                           
# Date: 2014-9-11                                                 
# Purpose: Extract station PBL data from ARL meteorological data   
# Note: Sample                                                    
#----------------------------------------------------------- 
from org.meteoinfo.data.meteodata import MeteoDataInfo 
import os.path 
from java.util import Date 
from java.util import Calendar 
from java.util import Locale 
from java.text import SimpleDateFormat 
 
#---- Set directories 
dataDir = 'P:/MetData/MetData/2011/' 
 
#---- Set times - this case show seasonal variation 
cal = Calendar.getInstance() 
cal.set(2011, 10, 1) 
sTime = cal.getTime() 
cal.set(2011, 11, 1) 
eTime = cal.getTime() 
aTime = sTime 
format = SimpleDateFormat('MMMyy', Locale.ENGLISH) 
format1 = SimpleDateFormat('yyyy-MM-dd HH:00') 
format2 = SimpleDateFormat('yyyyMMddHH') 
cal.setTime(sTime) 
 
#---- Create a MeteoDataInfo object 
mydata = MeteoDataInfo() 
#---- Set lon/lat of the location 
lon = 20.5 
lat = 40.5 
#---- Set output text file 
outpath = os.path.join(dataDir, 'PBL_data-1.txt') 
outf = open(outpath, 'w') 
outf.write('Time,PBLH') 
#---- Loop 
i = 0 
sum = 0.0 
while aTime <= eTime: 
          for w in range(5): 
                  inFile = dataDir + "gdas1." + format.format(aTime) + ".w" + str(w + 1) 
                print inFile 
                if os.path.isfile(inFile): 
                        mydata.openARLData(inFile)  #Open ARL data file 
                        tNum = mydata.getDataInfo().getTimeNum()    
                        print tNum 
                        for t in range(tNum): 
                                dTime = mydata.getDataInfo().getTimes().get(t) 
                                print format1.format(dTime) 
                                mydata.setTimeIndex(t)   #Set time index 
                                gData = mydata.getGridData("PBLH")  #Get PBLH grid data 
                                pblh = gData.toStation(lon, lat)  #Intepolate PBLH data to the location 
                                print 'PBLH = %10.2f' %pblh 
                                outf.write(os.linesep) 
                                outf.write(format2.format(dTime) + ',%-7.2f ' %pblh) 
                                sum += pblh 
                                i += 1 
                           
        cal.add(Calendar.MONTH, 1) 
        aTime = cal.getTime() 
 
outf.close() 
#---- Calculate seasonal average PBLH value 
ave = sum / i         
print 'Average PBLH: %10.2f' %ave 
 
print 'Finished!'

#--------------------------------------------------------
# Author: Yaqiang Wang
# Date: 2014-9-11
# Purpose: Extract station wind data from ARL meteorological data
# Note: Sample
#-----------------------------------------------------------
from org.meteoinfo.data import DataMath
from org.meteoinfo.data.meteodata import MeteoDataInfo
import os.path
from java.util import Date
from java.util import Calendar
from java.util import Locale
from java.text import SimpleDateFormat
#---- Set directories
dataDir = 'P:/MetData/MetData/2011/'
#---- Set times - this case show seasonal variation
cal = Calendar.getInstance()
cal.set(2011, 10, 1)
sTime = cal.getTime()
cal.set(2011, 10, 10)
eTime = cal.getTime()
aTime = sTime
format = SimpleDateFormat('MMMyy', Locale.ENGLISH)
format1 = SimpleDateFormat('yyyy-MM-dd HH:00')
format2 = SimpleDateFormat('yyyyMMddHH')
cal.setTime(sTime)
#---- Create a MeteoDataInfo object
mydata = MeteoDataInfo()
#---- Set lon/lat of the location
lon = 20.5
lat = 40.5
#---- Set output text file
outpath = os.path.join(dataDir, 'wind_data.txt')
outf = open(outpath, 'w')
outf.write('Time,Pressure,WindDir,WindSpeed')
#---- Loop
while aTime <= eTime:
for w in range(5):
inFile = dataDir + 'gdas1.' + format.format(aTime) + '.w' + str(w + 1)
print inFile
if os.path.isfile(inFile):
mydata.openARLData(inFile) #Open ARL data file
dataInfo = mydata.getDataInfo()
tNum = dataInfo.getTimeNum() #Get time number
print tNum
zdim = dataInfo.getVariable('UWND').getZDimension()
zNum = zdim.getDimLength() #Get level number
for t in range(tNum):
dTime = dataInfo.getTimes().get(t)
print format1.format(dTime)
mydata.setTimeIndex(t) #Set time index
for z in range(zNum):
mydata.setLevelIndex(z) #Set level index
lev = zdim.getDimValue().get(z)
uData = mydata.getGridData('UWND') #Get U grid data
vData = mydata.getGridData('VWND') #Get V grid data
#Calculate wind direction and wind speed from U/V
windData = DataMath.getDSFromUV(uData, vData)
wdData = windData[0]
wsData = windData[1]
wd = wdData.toStation(lon, lat) #Intepolate grid data to the location
ws = wsData.toStation(lon, lat)
print 'Wind direction = %10.2f; Wind speed = %10.2f' %(wd,ws)
outf.write(os.linesep)
line = format2.format(dTime) + ',%d'%lev + ',%.2f '%wd + ',%.2f'%ws
outf.write(line)
cal.add(Calendar.MONTH, 1)
aTime = cal.getTime()
outf.close()
print 'Finished!'

#--------------------------------------------------------         
# Author: Yaqiang Wang                                           
# Date: 2014-9-11                                                 
# Purpose: Extract station wind data from ARL meteorological data   
# Note: Sample                                                    
#----------------------------------------------------------- 
from org.meteoinfo.data import DataMath 
from org.meteoinfo.data.meteodata import MeteoDataInfo 
import os.path 
from java.util import Date 
from java.util import Calendar 
from java.util import Locale 
from java.text import SimpleDateFormat 
 
#---- Set directories 
dataDir = 'P:/MetData/MetData/2011/' 
 
#---- Set times - this case show seasonal variation 
cal = Calendar.getInstance() 
cal.set(2011, 10, 1) 
sTime = cal.getTime() 
cal.set(2011, 10, 10) 
eTime = cal.getTime() 
aTime = sTime 
format = SimpleDateFormat('MMMyy', Locale.ENGLISH) 
format1 = SimpleDateFormat('yyyy-MM-dd HH:00') 
format2 = SimpleDateFormat('yyyyMMddHH') 
cal.setTime(sTime) 
 
#---- Create a MeteoDataInfo object 
mydata = MeteoDataInfo() 
#---- Set lon/lat of the location 
lon = 20.5 
lat = 40.5 
#---- Set output text file 
outpath = os.path.join(dataDir, 'wind_data.txt') 
outf = open(outpath, 'w') 
outf.write('Time,Pressure,WindDir,WindSpeed') 
#---- Loop 
while aTime <= eTime: 
          for w in range(5): 
                  inFile = dataDir + 'gdas1.' + format.format(aTime) + '.w' + str(w + 1) 
                print inFile 
                if os.path.isfile(inFile): 
                        mydata.openARLData(inFile)  #Open ARL data file 
                        dataInfo = mydata.getDataInfo() 
                        tNum = dataInfo.getTimeNum()    #Get time number 
                        print tNum 
                        zdim = dataInfo.getVariable('UWND').getZDimension() 
                        zNum = zdim.getDimLength()    #Get level number 
                        for t in range(tNum): 
                                dTime = dataInfo.getTimes().get(t) 
                                print format1.format(dTime) 
                                mydata.setTimeIndex(t)   #Set time index 
                                for z in range(zNum): 
                                        mydata.setLevelIndex(z)    #Set level index 
                                        lev = zdim.getDimValue().get(z) 
                                        uData = mydata.getGridData('UWND')    #Get U grid data 
                                        vData = mydata.getGridData('VWND')    #Get V grid data 
                                        #Calculate wind direction and wind speed from U/V 
                                        windData = DataMath.getDSFromUV(uData, vData) 
                                        wdData = windData[0] 
                                        wsData = windData[1] 
                                        wd = wdData.toStation(lon, lat)  #Intepolate grid data to the location 
                                        ws = wsData.toStation(lon, lat) 
                                        print 'Wind direction = %10.2f; Wind speed = %10.2f' %(wd,ws) 
                                        outf.write(os.linesep) 
                                        line = format2.format(dTime) + ',%d'%lev + ',%.2f '%wd + ',%.2f'%ws 
                                        outf.write(line) 
                           
        cal.add(Calendar.MONTH, 1) 
        aTime = cal.getTime() 
 
outf.close() 
 
print 'Finished!'

qxtlyf · 发表于 2014-9-14 19:59:31

学习了，谢谢王老师

qx124 · 发表于 2014-9-14 20:37:58

好东西，O(∩_∩)O谢谢分享

finalkiller00 · 发表于 2014-9-14 22:26:24

{:eb303:}{:eb303:}

yrovl · 发表于 2014-9-15 08:57:54

学习了，谢谢王老师！

zhangyanlin · 发表于 2015-5-30 15:08:13

王老师，运行之后，出现错误代码，如下：
File "<iostream>", line 42
? ?? ?? ? for w in range(5):
^
SyntaxError: no viable alternative at character '?'

zhangyanlin · 发表于 2015-5-30 15:36:31

版主，运行后出现以下错误代码：
File "", line 42
          for w in range(5):

^
SyntaxError: expected an indented block

zhangyanlin · 发表于 2015-5-30 16:25:30

代码改了一下可以了，缩进之前有点问题。

#--------------------------------------------------------
# Author: Yaqiang Wang
# Date: 2014-9-11
# Purpose: Extract station wind data from ARL meteorological data
# Note: Sample
#-----------------------------------------------------------
from org.meteoinfo.data import DataMath
from org.meteoinfo.data.meteodata import MeteoDataInfo
import os.path
from java.util import Date
from java.util import Calendar
from java.util import Locale
from java.text import SimpleDateFormat
#---- Set directories
dataDir = 'G:\\Meteodata\\2012\\'
#---- Set times - this case show seasonal variation
cal = Calendar.getInstance()
cal.set(2012, 10, 1)
sTime = cal.getTime()
cal.set(2012, 10, 10)
eTime = cal.getTime()
aTime = sTime
format = SimpleDateFormat('MMMyy', Locale.ENGLISH)
format1 = SimpleDateFormat('yyyy-MM-dd HH:00')
format2 = SimpleDateFormat('yyyyMMddHH')
cal.setTime(sTime)
#---- Create a MeteoDataInfo object
mydata = MeteoDataInfo()
#---- Set lon/lat of the location
lon = 20.5
lat = 40.5
#---- Set output text file
outpath = os.path.join(dataDir, 'wind_data.txt')
outf = open(outpath, 'w')
outf.write('Time,Pressure,WindDir,WindSpeed')
#---- Loop
while aTime <= eTime:
for w in range(5):
inFile = dataDir + 'gdas1.' + format.format(aTime) + '.w' + str(w + 1)
print inFile
if os.path.isfile(inFile):
mydata.openARLData(inFile) #Open ARL data file
dataInfo = mydata.getDataInfo()
tNum = dataInfo.getTimeNum() #Get time number
print tNum
zdim = dataInfo.getVariable('UWND').getZDimension()
zNum = zdim.getDimLength() #Get level number
for t in range(tNum):
dTime = dataInfo.getTimes().get(t)
print format1.format(dTime)
mydata.setTimeIndex(t) #Set time index
for z in range(zNum):
mydata.setLevelIndex(z) #Set level index
lev = zdim.getDimValue().get(z)
uData = mydata.getGridData('UWND') #Get U grid data
vData = mydata.getGridData('VWND') #Get V grid data
#Calculate wind direction and wind speed from U/V
windData = DataMath.getDSFromUV(uData, vData)
wdData = windData[0]
wsData = windData[1]
wd = wdData.toStation(lon, lat) #Intepolate grid data to the location
ws = wsData.toStation(lon, lat)
print 'Wind direction = %10.2f; Wind speed = %10.2f' %(wd,ws)
outf.write(os.linesep)
line = format2.format(dTime) + ',%d'%lev + ',%.2f '%wd + ',%.2f'%ws
outf.write(line)
cal.add(Calendar.MONTH, 1)
aTime = cal.getTime()
outf.close()
print 'Finished!'

复制代码

精卫 · 发表于 2015-5-30 20:20:48

这是个很实用的东西啊王老师我才看见，用过由格点转站点的，记得是用surfer完成的，看到您的这个东西恒感慨，虽然也会写点程序什么的，但是对您的脚本确是一窍不通，我是您软件最早的试用者，却一直不会您的脚本，汗颜啊

yagudu · 发表于 2016-7-19 20:15:08

王老师，求边界层高度的我用了一下，出现一个错误，不知道是怎么回事
>>> run script...
File "<iostream>", line 40
print inFile
^
IndentationError: unindent does not match any outer indentation level

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

MeteoInfo脚本示例 - 从再分析数据中提取站点数据

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

评分