MeteoInfo脚本教程（十七）：站点数据分析

MeteoInfo · 发表于 2014-12-21 22:30:36

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

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

x

站点数据也有类似格点数据的四则运算等，用法也很类似。此外站点数据经常需要插值为格点数据然后进行等值线等分析，因此插值算法也很重要。MeteoInfo站点数据插值为格点数据的算法主要有两种：反距离权法（IDW）和Cressman方法。其中IDW方法又分为IDW_Radius和IDW_Neighbour两种。IDW_Radius和Cressman插值方法都需要搜索半径，在站点稀疏的地区可能会插出缺测值（搜索半径内没有站点）。设置插值方法和参数的类是InterpolationSetting类，位于org.meteoinfo.geoprocess.analysis包中，其构造函数有两个：

IDW方法：
InterpolationSetting(double minX, double maxX, double minY, double maxY, int xNum, int yNum, String aInterMethod, float radius, int minNum)

Cressman方法：
InterpolationSetting(double minX, double maxX, double minY, double maxY, int xNum, int yNum, String aInterMethod, List<Double> radList)

前六个参数定义了格点场的范围和分辨率，第7个参数是方法名（'IDW_Radius', 'IDW_Neighbour', 'Cressman'），radius是搜索半径，minNum是使用的最少站点数，Cressman插值方法需要多个搜索半径（radList参数）。

这里给出一个计算站点累计降水并插值成格点绘制等值线填色图的例子。内容比较丰富，包括了站点数据运算（相加）、插值、添加南海脚图、投影、标注等功能，希望大家能仔细理解。

# encoding=utf-8
#-----------------------------------------------------
# Author: Yaqiang Wang
# Date: 2014-12-6
# Purpose: Interpolate station data to grid data
# Note: Sample
#-----------------------------------------------------
#---- Import classes
print 'Import classes...'
from org.meteoinfo.layout import MapLayout
from org.meteoinfo.data import StationData
from org.meteoinfo.data.mapdata import MapDataManage
from org.meteoinfo.data.meteodata import MeteoDataInfo, DrawMeteoData
from org.meteoinfo.legend import LegendScheme, LegendType, PolygonBreak, MapFrame
from org.meteoinfo.shape import ShapeTypes
from org.meteoinfo.layout import LegendStyles, LayoutMap
from org.meteoinfo.global import Extent
from org.meteoinfo.projection import ProjectionInfo
from org.meteoinfo.geoprocess.analysis import InterpolationSetting
import os
from datetime import datetime
from datetime import timedelta
from java.awt import Color
from java.awt import Font
from javax.swing import JFrame
#---- Set directories
print 'Set directories...'
baseDir = 'D:/MyProgram/Distribution/java/MeteoInfo/MeteoInfo'
dataDir = os.path.join(baseDir, 'sample/MICAPS')
mapDir = os.path.join(baseDir, 'map')
#---- Create MapLayout object
mapLayout = MapLayout()
mapFrame = mapLayout.getActiveMapFrame()
mapView = mapFrame.getMapView()
layoutMap = mapLayout.getActiveLayoutMap()
#---- Load layers
print 'Load layers...'
bou2Layer = MapDataManage.loadLayer(os.path.join(mapDir, 'bou2_4p.shp'))
lb = bou2Layer.getLegendScheme().getLegendBreaks().get(0)
lb.setDrawFill(False)
lb.setOutlineColor(Color.gray)
mapFrame.addLayer(bou2Layer)
bou1Layer = MapDataManage.loadLayer(os.path.join(mapDir, 'bou1_4l.shp'))
lb = bou1Layer.getLegendScheme().getLegendBreaks().get(0)
lb.setColor(Color.blue)
mapFrame.addLayer(bou1Layer)
bou1Layer_1 = MapDataManage.loadLayer(os.path.join(mapDir, 'bou1_4l.shp'))
lb = bou1Layer_1.getLegendScheme().getLegendBreaks().get(0)
lb.setColor(Color.blue)
chinaLayer = MapDataManage.loadLayer(os.path.join(mapDir, 'china.shp'))
chinaLayer.setVisible(False)
mapFrame.addLayer(chinaLayer)
cityLayer = MapDataManage.loadLayer(os.path.join(mapDir, 'res1_4m.shp'))
lb = cityLayer.getLegendScheme().getLegendBreaks().get(0)
lb.setSize(6)
lb.setColor(Color.red)
labset = cityLayer.getLabelSet()
labset.setFieldName('NAME')
font = Font(u'楷体', Font.PLAIN, 16)
labset.setLabelFont(font)
labset.setYOffset(15)
cityLayer.addLabels()
mapFrame.addLayer(cityLayer)
#---- Add South China Sea
aMapFrame = MapFrame()
aLayoutMap = LayoutMap(aMapFrame)
aLayoutMap.setDrawGridLabel(False)
aLayoutMap.setDrawGridTickLine(False)
aLayoutMap.setLeft(40)
aLayoutMap.setTop(350)
aLayoutMap.setWidth(85)
aLayoutMap.setHeight(109)
mapLayout.addElement(aLayoutMap)
aMapFrame.addLayer(bou1Layer_1)
aProjInfo = ProjectionInfo("+proj=lcc+lat_1=25+lat_2=47+lon_0=105")
aMapFrame.getMapView().projectLayers(aProjInfo)
aMapFrame.getMapView().zoomToExtentLonLatEx(Extent(106.5,122.5,1,23))
#---- Project mapview
print 'Project mapview...'
projStr = '+proj=lcc+lat_1=25+lat_2=47+lon_0=105'
projInfo = ProjectionInfo(projStr)
mapView.projectLayers(projInfo)
extent = Extent(78,130,15,53)
mapView.zoomToExtentLonLatEx(extent)
#---- Create MeteoDataInfo object
mdi = MeteoDataInfo()
#---- Set start/end time
stime = datetime(2010,10,14,14)
etime = datetime(2010,10,14,20)
#---- Loop
print 'Get sum station data...'
sdata = StationData()
atime = stime
i = 0
while atime < etime:
#---- Open a MICAPS data file
fn = os.path.join(dataDir, stime.strftime('%y%m%d%H') + '.000')
mdi.openMICAPSData(fn)
#---- Sum precipitation station data
if i == 0:
sdata = mdi.getStationData('Precipitation6h')
else:
sdata = sdata.add(mdi.getStationData('Precipitation6h'))
atime = atime + timedelta(hours=6)
#---- Interpolate station data to grid data
print 'Interpolate station data to grid data...'
interpSet = InterpolationSetting(60,140,-20,60,160,160,"IDW_Radius",2,1)
#radList = [10.0, 8.0, 6.0, 4.0]
#interpSet = InterpolationSetting(60,140,-20,60,160,160,"Cressman",radList)
gdata = sdata.interpolateData(interpSet)
#---- Set legend scheme
ls = LegendScheme(ShapeTypes.Polygon)
ls.setLegendType(LegendType.GraduatedColor)
values = [0,0.1,1,2,5,10,20,25,50,100,1000]
colors = [Color(255,255,255),Color(170,240,255),Color(120,230,240),Color(200,220,50),Color(240,220,20),
Color(255,120,10),Color(255,90,10),Color(240,40,0),Color(180,10,0),Color(120,10,0)]
lbs = ls.getLegendBreaks()
for i in range(0, len(colors)):
lb = PolygonBreak()
lb.setColor(colors)
lb.setStartValue(values)
lb.setEndValue(values[i + 1])
if i == 0:
lb.setCaption('< ' + str(values[i + 1]))
elif i == len(colors) - 1:
lb.setCaption('> ' + str(values))
else:
lb.setCaption(str(values) + ' - ' + str(values[i + 1]))
lb.setDrawOutline(False)
lbs.add(lb)
#---- Create shaded layer
print 'Create shaded layer...'
layer = DrawMeteoData.createShadedLayer(gdata, ls, 'Rain_shaded', 'Rain', True)
layer.setMaskout(True)
#---- Add layer
mapFrame.addLayer(layer)
mapFrame.moveLayer(layer, 0)
#---- Add title
title = mapLayout.addText(u'全国降水量实况图', 280, 40, u'黑体', 18)
stime = stime + timedelta(hours=-6)
subTitle = mapLayout.addText(u'（' + stime.strftime('%Y-%m-%d %H:00') +
u' 至 ' + etime.strftime('%Y-%m-%d %H:00') + u'）', 280, 60, u'黑体', 16)
#---- Set layout map
print 'Set layout map...'
layoutMap.setDrawGridLine(False)
layoutMap.setDrawNeatLine(False)
layoutMap.setDrawGridLabel(False)
layoutMap.setDrawGridTickLine(False)
layoutMap.setLeft(10)
layoutMap.setTop(10)
layoutMap.setWidth(620)
layoutMap.setHeight(450)
#---- Set maskout
mapView.getMaskOut().setMask(True)
mapView.getMaskOut().setMaskLayer(chinaLayer.getLayerName())
#---- Set mapframe
mapFrame.setGridXDelt(10)
mapFrame.setGridYDelt(10)
#---- Add legend
legend = mapLayout.addLegend(575, 250)
legend.setLegendStyle(LegendStyles.Normal)
legend.setLegendLayer(layer)
legend.setFont(Font(u'宋体', Font.PLAIN, 12))
legend.setTitle(u'降水量（毫米）')
frame = JFrame('MeteoInfo Script Sample', size = (750, 530))
frame.add(mapLayout)
frame.visible = True
print 'Finished!'

# encoding=utf-8 
#----------------------------------------------------- 
# Author: Yaqiang Wang 
# Date: 2014-12-6 
# Purpose: Interpolate station data to grid data 
# Note: Sample 
#----------------------------------------------------- 
#---- Import classes 
print 'Import classes...' 
from org.meteoinfo.layout import MapLayout 
from org.meteoinfo.data import StationData 
from org.meteoinfo.data.mapdata import MapDataManage 
from org.meteoinfo.data.meteodata import MeteoDataInfo, DrawMeteoData 
from org.meteoinfo.legend import LegendScheme, LegendType, PolygonBreak, MapFrame 
from org.meteoinfo.shape import ShapeTypes 
from org.meteoinfo.layout import LegendStyles, LayoutMap 
from org.meteoinfo.global import Extent 
from org.meteoinfo.projection import ProjectionInfo 
from org.meteoinfo.geoprocess.analysis import InterpolationSetting 
import os 
from datetime import datetime 
from datetime import timedelta 
from java.awt import Color 
from java.awt import Font 
from javax.swing import JFrame 
 
#---- Set directories 
print 'Set directories...' 
baseDir = 'D:/MyProgram/Distribution/java/MeteoInfo/MeteoInfo' 
dataDir = os.path.join(baseDir, 'sample/MICAPS') 
mapDir = os.path.join(baseDir, 'map') 
 
#---- Create MapLayout object 
mapLayout = MapLayout() 
mapFrame = mapLayout.getActiveMapFrame() 
mapView = mapFrame.getMapView() 
layoutMap = mapLayout.getActiveLayoutMap() 
 
#---- Load layers 
print 'Load layers...' 
bou2Layer = MapDataManage.loadLayer(os.path.join(mapDir, 'bou2_4p.shp')) 
lb = bou2Layer.getLegendScheme().getLegendBreaks().get(0) 
lb.setDrawFill(False) 
lb.setOutlineColor(Color.gray) 
mapFrame.addLayer(bou2Layer) 
bou1Layer = MapDataManage.loadLayer(os.path.join(mapDir, 'bou1_4l.shp')) 
lb = bou1Layer.getLegendScheme().getLegendBreaks().get(0) 
lb.setColor(Color.blue) 
mapFrame.addLayer(bou1Layer) 
bou1Layer_1 = MapDataManage.loadLayer(os.path.join(mapDir, 'bou1_4l.shp')) 
lb = bou1Layer_1.getLegendScheme().getLegendBreaks().get(0) 
lb.setColor(Color.blue) 
chinaLayer = MapDataManage.loadLayer(os.path.join(mapDir, 'china.shp')) 
chinaLayer.setVisible(False) 
mapFrame.addLayer(chinaLayer) 
cityLayer = MapDataManage.loadLayer(os.path.join(mapDir, 'res1_4m.shp')) 
lb = cityLayer.getLegendScheme().getLegendBreaks().get(0) 
lb.setSize(6) 
lb.setColor(Color.red) 
labset = cityLayer.getLabelSet() 
labset.setFieldName('NAME') 
font = Font(u'楷体', Font.PLAIN, 16) 
labset.setLabelFont(font) 
labset.setYOffset(15) 
cityLayer.addLabels() 
mapFrame.addLayer(cityLayer) 
 
#---- Add South China Sea 
aMapFrame = MapFrame() 
aLayoutMap = LayoutMap(aMapFrame) 
aLayoutMap.setDrawGridLabel(False) 
aLayoutMap.setDrawGridTickLine(False) 
aLayoutMap.setLeft(40) 
aLayoutMap.setTop(350) 
aLayoutMap.setWidth(85) 
aLayoutMap.setHeight(109) 
mapLayout.addElement(aLayoutMap) 
aMapFrame.addLayer(bou1Layer_1) 
aProjInfo = ProjectionInfo("+proj=lcc+lat_1=25+lat_2=47+lon_0=105") 
aMapFrame.getMapView().projectLayers(aProjInfo) 
aMapFrame.getMapView().zoomToExtentLonLatEx(Extent(106.5,122.5,1,23)) 
 
#---- Project mapview 
print 'Project mapview...' 
projStr = '+proj=lcc+lat_1=25+lat_2=47+lon_0=105' 
projInfo = ProjectionInfo(projStr) 
mapView.projectLayers(projInfo) 
extent = Extent(78,130,15,53) 
mapView.zoomToExtentLonLatEx(extent) 
 
#---- Create MeteoDataInfo object 
mdi = MeteoDataInfo() 
 
#---- Set start/end time 
stime = datetime(2010,10,14,14) 
etime = datetime(2010,10,14,20) 
 
#---- Loop 
print 'Get sum station data...' 
sdata = StationData() 
atime = stime 
i = 0 
while atime < etime: 
  #---- Open a MICAPS data file 
  fn = os.path.join(dataDir, stime.strftime('%y%m%d%H') + '.000') 
  mdi.openMICAPSData(fn) 
 
  #---- Sum precipitation station data 
  if i == 0: 
    sdata = mdi.getStationData('Precipitation6h') 
  else: 
    sdata = sdata.add(mdi.getStationData('Precipitation6h')) 
 
  atime = atime + timedelta(hours=6) 
   
#---- Interpolate station data to grid data 
print 'Interpolate station data to grid data...' 
interpSet = InterpolationSetting(60,140,-20,60,160,160,"IDW_Radius",2,1) 
#radList = [10.0, 8.0, 6.0, 4.0] 
#interpSet = InterpolationSetting(60,140,-20,60,160,160,"Cressman",radList) 
gdata = sdata.interpolateData(interpSet) 
 
#---- Set legend scheme 
ls = LegendScheme(ShapeTypes.Polygon) 
ls.setLegendType(LegendType.GraduatedColor) 
values = [0,0.1,1,2,5,10,20,25,50,100,1000] 
colors = [Color(255,255,255),Color(170,240,255),Color(120,230,240),Color(200,220,50),Color(240,220,20), 
  Color(255,120,10),Color(255,90,10),Color(240,40,0),Color(180,10,0),Color(120,10,0)] 
lbs = ls.getLegendBreaks() 
for i in range(0, len(colors)): 
  lb = PolygonBreak() 
  lb.setColor(colors) 
  lb.setStartValue(values) 
  lb.setEndValue(values[i + 1]) 
  if i == 0: 
    lb.setCaption('< ' + str(values[i + 1])) 
  elif i == len(colors) - 1: 
    lb.setCaption('> ' + str(values)) 
  else: 
    lb.setCaption(str(values) + ' - ' + str(values[i + 1])) 
  lb.setDrawOutline(False) 
  lbs.add(lb) 
 
#---- Create shaded layer 
print 'Create shaded layer...'   
layer = DrawMeteoData.createShadedLayer(gdata, ls, 'Rain_shaded', 'Rain', True) 
layer.setMaskout(True) 
 
#---- Add layer 
mapFrame.addLayer(layer) 
mapFrame.moveLayer(layer, 0) 
 
#---- Add title 
title = mapLayout.addText(u'全国降水量实况图', 280, 40, u'黑体', 18) 
stime = stime + timedelta(hours=-6) 
subTitle = mapLayout.addText(u'（' + stime.strftime('%Y-%m-%d %H:00') + 
  u' 至 ' + etime.strftime('%Y-%m-%d %H:00') + u'）', 280, 60, u'黑体', 16) 
 
#---- Set layout map 
print 'Set layout map...' 
layoutMap.setDrawGridLine(False) 
layoutMap.setDrawNeatLine(False) 
layoutMap.setDrawGridLabel(False) 
layoutMap.setDrawGridTickLine(False) 
layoutMap.setLeft(10) 
layoutMap.setTop(10) 
layoutMap.setWidth(620) 
layoutMap.setHeight(450) 
 
#---- Set maskout 
mapView.getMaskOut().setMask(True) 
mapView.getMaskOut().setMaskLayer(chinaLayer.getLayerName()) 
 
#---- Set mapframe 
mapFrame.setGridXDelt(10) 
mapFrame.setGridYDelt(10) 
 
#---- Add legend 
legend = mapLayout.addLegend(575, 250) 
legend.setLegendStyle(LegendStyles.Normal) 
legend.setLegendLayer(layer) 
legend.setFont(Font(u'宋体', Font.PLAIN, 12)) 
legend.setTitle(u'降水量（毫米）') 
 
frame = JFrame('MeteoInfo Script Sample', size = (750, 530)) 
frame.add(mapLayout) 
frame.visible = True 
print 'Finished!'

BDBD168 · 发表于 2014-12-22 08:38:59

不错，

{:eb502:}{:eb502:}{:eb502:}{:eb502:}

无敌泡 · 发表于 2014-12-24 09:37:25

楼主好厉害

嘴角上扬 · 发表于 2014-12-31 09:00:16

图画的好漂亮

yingkr · 发表于 2015-6-7 17:41:58

太厉害，介绍的很详细~~~~~

iij · 发表于 2015-12-22 16:00:03

想请问一下地图根据范经纬度围缩放在C#里面怎么使用呢？

wm358708635 · 发表于 2015-12-24 21:34:03

很好一啊哈哈哈

zfvv · 发表于 2015-12-30 12:55:36

不知道王老师在站点数据较多，需要绘制等值线时，是否对三角剖分追踪等值线算法实现？

MeteoInfo · 发表于 2015-12-30 16:06:07

zfvv 发表于 2015-12-30 12:55
不知道王老师在站点数据较多，需要绘制等值线时，是否对三角剖分追踪等值线算法实现？

没研究过那个算法。

kcuas · 发表于 2017-8-18 10:53:46

您好，我是用Java二次开发的，生成的图例的名称我按照

legend.setTitle("降水量（mm）");

复制代码

结果却显示contourlayer,这是为什么呢，要怎么解决

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

MeteoInfo脚本教程（十七）：站点数据分析

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

本帖被以下淘专辑推荐: