MeteoInfoLab脚本示例：站点数据读取插值

MeteoInfo

数据来自此贴：http://bbs.06climate.com/forum.p ... &extra=page%3D2

这里示例用MeteoInfoLab脚本读取站点数据并匹配经纬度，绘制散点图，插值为格点数据并绘制等值线图的过程。

1. #Set file name
   
2. stfn = 'E:/Temp/station.txt'
   
3. datafn = 'E:/Temp/201511_auto_hr.txt'
   
4. tstr = '2015110101'
   
5. #Read station id, lon, lat
   
6. st_table = readtable(stfn, headerlines=-1, format='%s%2f')
   
7. stids = st_table['Col_0']
   
8. lons = st_table['Col_1']
   
9. lats = st_table['Col_2']
   
10. #Read station temperature data at a specific time
    
11. stations = []
    
12. temp = []
    
13. lon = []
    
14. lat = []
    
15. dataf = open(datafn)
    
16. for line in dataf:
    
17.     a = line.split()
    
18.     if a[1] == tstr:
    
19.         if a[0] in stids:
    
20.             idx = stids.index(a[0])
    
21.             stations.append(a[0])
    
22.             temp.append(float(a[3]))
    
23.             lon.append(lons[idx])
    
24.             lat.append(lats[idx])
    
25. dataf.close()
    
26. 
    
27. temp = array(temp)
    
28. temp[temp==-9997] = nan
    
29. #Plot
    
30. axesm()
    
31. ltaiwan = shaperead('D:/Temp/Map/taiwan.shp')
    
32. geoshow(ltaiwan, edgecolor='k')
    
33. layer = scatterm(lon, lat, temp, 10)
    
34. colorbar(layer)
    
35. xlim(119.8, 122.2)
    
36. ylim(21.8, 25.5)

#Set file name<br /> stfn = 'E:/Temp/station.txt'<br /> datafn = 'E:/Temp/201511_auto_hr.txt'<br /> tstr = '2015110101'<br /> #Read station id, lon, lat<br /> st_table = readtable(stfn, headerlines=-1, format='%s%2f')<br /> stids = st_table['Col_0']<br /> lons = st_table['Col_1']<br /> lats = st_table['Col_2']<br /> #Read station temperature data at a specific time<br /> stations = []<br /> temp = []<br /> lon = []<br /> lat = []<br /> dataf = open(datafn)<br /> for line in dataf:<br />     a = line.split()<br />     if a[1] == tstr:<br />         if a[0] in stids:<br />             idx = stids.index(a[0])<br />             stations.append(a[0])<br />             temp.append(float(a[3]))<br />             lon.append(lons[idx])<br />             lat.append(lats[idx])<br /> dataf.close()<br /> <br /> temp = array(temp)<br /> temp[temp==-9997] = nan<br /> #Plot<br /> axesm()<br /> ltaiwan = shaperead('D:/Temp/Map/taiwan.shp')<br /> geoshow(ltaiwan, edgecolor='k')<br /> layer = scatterm(lon, lat, temp, 10)<br /> colorbar(layer)<br /> xlim(119.8, 122.2)<br /> ylim(21.8, 25.5)

1. #Set file name
   
2. stfn = 'E:/Temp/station.txt'
   
3. datafn = 'E:/Temp/201511_auto_hr.txt'
   
4. tstr = '2015110101'
   
5. #Read station id, lon, lat
   
6. st_table = readtable(stfn, headerlines=-1, format='%s%2f')
   
7. stids = st_table['Col_0']
   
8. lons = st_table['Col_1']
   
9. lats = st_table['Col_2']
   
10. #Read station temperature data at a specific time
    
11. stations = []
    
12. temp = []
    
13. lon = []
    
14. lat = []
    
15. dataf = open(datafn)
    
16. for line in dataf:
    
17.     a = line.split()
    
18.     if a[1] == tstr:
    
19.         if a[0] in stids:
    
20.             idx = stids.index(a[0])
    
21.             stations.append(a[0])
    
22.             temp.append(float(a[3]))
    
23.             lon.append(lons[idx])
    
24.             lat.append(lats[idx])
    
25. dataf.close()
    
26. 
    
27. temp = array(temp)
    
28. temp[temp==-9997] = nan
    
29. #Interpolate to grid data
    
30. lon = array(lon)
    
31. lat = array(lat)
    
32. x = arange(119.8, 122.2, 0.1)
    
33. y = arange(21.8, 25.5, 0.1)
    
34. gtemp,xx,yy = griddata([lon, lat], temp, xi=[x, y], method='idw')
    
35. #Plot
    
36. axesm()
    
37. ltaiwan = shaperead('D:/Temp/Map/taiwan.shp')
    
38. tw_city = shaperead('D:/Temp/Map/taiwan_city.shp')
    
39. tw_xc = shaperead('D:/Temp/Map/taiwan_xc.shp')
    
40. geoshow(ltaiwan, edgecolor='gray')
    
41. geoshow(tw_xc, size=4, edgecolor='k', labelfield='NAME', \
    
42.     fontname=u'楷体', fontsize=14, yoffset=15)
    
43. geoshow(tw_city, facecolor='b', edgecolor='k', size=6, labelfield='NAME', \
    
44.     fontname=u'楷体', fontsize=16, yoffset=15)
    
45. tw_city.movelabel(u'台北', 0, -25)
    
46. layer = contourfm(xx, yy, gtemp, 20)
    
47. masklayer(ltaiwan, [layer])
    
48. colorbar(layer)
    
49. xlim(119.8, 122.2)
    
50. ylim(21.8, 25.5)

#Set file name<br /> stfn = 'E:/Temp/station.txt'<br /> datafn = 'E:/Temp/201511_auto_hr.txt'<br /> tstr = '2015110101'<br /> #Read station id, lon, lat<br /> st_table = readtable(stfn, headerlines=-1, format='%s%2f')<br /> stids = st_table['Col_0']<br /> lons = st_table['Col_1']<br /> lats = st_table['Col_2']<br /> #Read station temperature data at a specific time<br /> stations = []<br /> temp = []<br /> lon = []<br /> lat = []<br /> dataf = open(datafn)<br /> for line in dataf:<br />     a = line.split()<br />     if a[1] == tstr:<br />         if a[0] in stids:<br />             idx = stids.index(a[0])<br />             stations.append(a[0])<br />             temp.append(float(a[3]))<br />             lon.append(lons[idx])<br />             lat.append(lats[idx])<br /> dataf.close()<br /> <br /> temp = array(temp)<br /> temp[temp==-9997] = nan<br /> #Interpolate to grid data <br /> lon = array(lon)<br /> lat = array(lat)<br /> x = arange(119.8, 122.2, 0.1)<br /> y = arange(21.8, 25.5, 0.1)<br /> gtemp,xx,yy = griddata([lon, lat], temp, xi=[x, y], method='idw') <br /> #Plot<br /> axesm()<br /> ltaiwan = shaperead('D:/Temp/Map/taiwan.shp')<br /> tw_city = shaperead('D:/Temp/Map/taiwan_city.shp')<br /> tw_xc = shaperead('D:/Temp/Map/taiwan_xc.shp')<br /> geoshow(ltaiwan, edgecolor='gray')<br /> geoshow(tw_xc, size=4, edgecolor='k', labelfield='NAME', \<br />     fontname=u'楷体', fontsize=14, yoffset=15)<br /> geoshow(tw_city, facecolor='b', edgecolor='k', size=6, labelfield='NAME', \<br />     fontname=u'楷体', fontsize=16, yoffset=15)<br /> tw_city.movelabel(u'台北', 0, -25)<br /> layer = contourfm(xx, yy, gtemp, 20)<br /> masklayer(ltaiwan, [layer])<br /> colorbar(layer)<br /> xlim(119.8, 122.2)<br /> ylim(21.8, 25.5)

随缘

{:eb502:}{:eb502:}

addname

哇，看起來好漂亮喔！
可是我沒有MeteoInfo這個程式。

addname

哇，看起來好漂亮喔！
可是我沒有MeteoInfo這個程式。

MeteoInfo

addname 发表于 2016-3-23 10:43
哇，看起來好漂亮喔！
可是我沒有MeteoInfo這個程式。

可以在这里免费下载使用：http://www.meteothinker.com/

addname

MeteoInfo 发表于 2016-3-23 10:52
可以在这里免费下载使用：http://www.meteothinker.com/

好的，謝謝您，來下載用看看！

clare

腻害！！！请教王老师，taiwan_xc.shp文件是怎么制作的？

MeteoInfo

clare 发表于 2016-3-24 11:06
腻害！！！请教王老师，taiwan_xc.shp文件是怎么制作的？

从全国县城shape文件选择性输出的

要吃早饭

{:eb502:}

想个名字

{:eb502:}{:eb502:}正在学习画散点图