MeteoInfoLab脚本示例：站点数据绘制等值线

MeteoInfo

gjj123 发表于 2016-9-13 09:53
老师您好，我编程的目标是想实现读取文本格式的数据，实现插值，于是我对您上面的例子进行了变动这样做了
...

没看明白到底想问什么？

gjj123

MeteoInfo 发表于 2016-9-13 10:00
没看明白到底想问什么？

主要是想实现读取文本文件，文件格式如下，包含精度，纬度，和值
#Set data folders
# coding=utf-8
basedir = 'C:/Users/jingjing_gu/Desktop/MeteoInfo_Java_1.3.4R3_Files/MeteoInfo'
datadir = os.path.join(basedir, 'sample/MICAPS')
mapdir = os.path.join(basedir, 'map')
#Read shape files
bou2_layer = shaperead(os.path.join(mapdir, 'bou2_4p.shp'))
bou1_layer = shaperead(os.path.join(mapdir, 'bou1_4l.shp'))
china_layer = shaperead(os.path.join(mapdir, 'china.shp'))
city_layer = shaperead(os.path.join(mapdir, 'res1_4m.shp'))
#Read station data

fn = 'E:/Temp/txt.dat'
ncol = numasciicol(fn)
nrow = numasciirow(fn)
a = asciiread(fn,shape=(nrow,ncol))
lon = a[:,0]
lat = a[:,1]
v = a[:,2]
axesm()

#griddata function - interpolate
x = arange(75, 135, 0.5)
y = arange(18, 55, 0.5)
prg = griddata((lon, lat), v, xi=(x, y), method='idw', convexhull=True)[0]

#Plot
axesm()
geoshow(bou2_layer, edgecolor='lightgray')
geoshow(bou1_layer, facecolor=(0,0,255))
geoshow(city_layer, facecolor='r', size=4, labelfield='NAME', fontname=u'楷体', fontsize=16, yoffset=15)
geoshow(china_layer, visible=False)
levs = [0.1, 1, 2, 5, 10, 20, 25, 50, 100]
cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \
    (240,40,0),(180,10,0),(120,10,0)]
layer = contourfm(prg, levs, colors=cols)
masklayer(china_layer, [layer])
colorbar(layer)
xlim(72, 136)
ylim(16, 55)
text(95, 52, u'全国PM2.5实况图', fontname=u'黑体', fontsize=16)
text(95, 50, u'(2010-10-14 08:00 至 2010-10-14 14:00)', fontname=u'黑体', fontsize=14)
#Add south China Sea
sc_layer = bou1_layer.clone()
axesm(position=[0.14,0.18,0.15,0.2], axison=False)
geoshow(sc_layer, facecolor=(0,0,255))
xlim(106, 123)
ylim(2, 23)
savefig("E:/Temp/ABC.BMP")
可是出来的结果不对，显示的点数并不是我txt.bat中所设的点

MeteoInfo

gjj123 发表于 2016-9-13 10:07
主要是想实现读取文本文件，文件格式如下，包含精度，纬度，和值
#Set data folders
# coding=utf-8

有错误信息吗？

gjj123

MeteoInfo 发表于 2016-9-13 10:11
有错误信息吗？

没有，运行时正确的，就是没有插值结果C:\Users\jingjing_gu\Desktop\捕获.JPG

MeteoInfo

gjj123 发表于 2016-9-13 10:17
没有，运行时正确的，就是没有插值结果

看不到你的图片

gjj123

MeteoInfo 发表于 2016-9-13 10:19
看不到你的图片

原始数据
经度       纬度     值
111.36   29.58   60
120.299 31.568  39
120.86   32.01   45
98         40        98
100       50        69
120.45   32.57   98
119.82  32.93    63
119.16  34.59    40
118.83  32.36    80
117.2    34.26    78


程序:
#Set data folders
# coding=utf-8
basedir = 'C:/Users/jingjing_gu/Desktop/MeteoInfo_Java_1.3.4R3_Files/MeteoInfo'
datadir = os.path.join(basedir, 'sample/MICAPS')
mapdir = os.path.join(basedir, 'map')
#Read shape files
bou2_layer = shaperead(os.path.join(mapdir, 'bou2_4p.shp'))
bou1_layer = shaperead(os.path.join(mapdir, 'bou1_4l.shp'))
china_layer = shaperead(os.path.join(mapdir, 'china.shp'))
city_layer = shaperead(os.path.join(mapdir, 'res1_4m.shp'))
#Read station data
fn = 'E:/Temp/txt.dat'
ncol = numasciicol(fn)
nrow = numasciirow(fn)
a = asciiread(fn,shape=(nrow,ncol))
lon = a[:,0]
lat = a[:,1]
v = a[:,2]
axesm()
#griddata function - interpolate
x = arange(75, 135,0.5)
y = arange(18, 55, 0.5)
prg = griddata((lon, lat), v, xi=(x, y),  method='idw', convexhull=True)[0]
#Plot
geoshow(bou2_layer, edgecolor='lightgray')
geoshow(bou1_layer, facecolor=(0,0,255))
geoshow(city_layer, facecolor='r', size=4, labelfield='NAME', fontname=u'楷体', fontsize=16, yoffset=15)
geoshow(china_layer, visible=False)
levs = [0.1, 1, 2, 5, 10, 20, 25, 50, 100]
cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \
    (240,40,0),(180,10,0),(120,10,0)]
layer = contourfm(prg, levs, colors=cols)
masklayer(china_layer, [layer])
colorbar(layer)
xlim(72, 136)
ylim(16, 55)
text(95, 52, u'全国PM2.5实况图', fontname=u'黑体', fontsize=16)
text(95, 50, u'(2010-10-14 08:00 至 2010-10-14 14:00)', fontname=u'黑体', fontsize=14)
#Add south China Sea
sc_layer = bou1_layer.clone()
axesm(position=[0.14,0.18,0.15,0.2], axison=False)
geoshow(sc_layer, facecolor=(0,0,255))
xlim(106, 123)
ylim(2, 23)
出的插值图感觉不合理，是不是程序中参数设置的有问题，请老师指教。

MeteoInfo

gjj123 发表于 2016-9-13 14:03
原始数据
经度       纬度     值
111.36   29.58   60

站点太少，分布又不均匀，不建议进行插值

gjj123

MeteoInfo 发表于 2016-9-13 14:07
站点太少，分布又不均匀，不建议进行插值

，谢谢！

gjj123

老师您好！关于插值还是有一些不明白的地方，还请多多指教
图1为采用成都的PM2.5的数据做的插值，插值结果并不理想，通过改变插值半径效果也不是很好

图1 通过设定半径参数插值结果图

但是我用ArcGis做出来的插值结果就比较好，ArcGis插值参数如图2所示，插值结果如图3所示。

图2

图3

是不是插值半径有问题，还是其他参数设置，ArcGis的插值半径通过指定点数是可变的？

MeteoInfo

gjj123 发表于 2016-9-18 09:19
老师您好！关于插值还是有一些不明白的地方，还请多多指教
图1为采用成都的PM2.5的数据做的插值 ...

把你用的数据文件和脚本程序附上，抽空帮你看看。