MeteoInfoLab脚本示例：线性拟合

MeteoInfo

MeteoInfoLab提供一个线性拟合函数linregress，参数是参与拟合的两个数据序列，返回拟合的斜率、截距和相关系数。有了上述拟合参数可以用polyval函数生成拟合数据（直线）。然后可以将数据、拟合线、公式等绘图。



脚本程序：

1. fn = os.path.join('D:/KeyData/PMMUL/data/54500_PMMUL_DA.csv')
   
2. if os.path.exists(fn):
   
3.     print fn
   
4.     tdata = readtable(fn, delimiter=',', format='%{yyyyMMdd}D%f%f%f%i%i%i')        
   
5.     pm10 = tdata['PM10']
   
6.     pm25 = tdata['PM2.5']
   
7.     slope, intercept, r, p, stderr, n = linregress(pm10, pm25)
   
8.     x = arange(0, 800, 100)
   
9.     y = polyval([slope, intercept], x)
   
10. 
    
11.     scatter(pm10, pm25, s=2, color='k', label='PM')
    
12.     plot(x, y, 'r-', linewidth=2)
    
13.     xlim(0, 800)
    
14.     ylim(0, 800)
    
15.     xlabel(r'$\rm{PM_{10}} \ (\mu g \ m^{-3})$')
    
16.     ylabel(r'$\rm{PM_{2.5}} \ (\mu g \ m^{-3})$')
    
17.     text(100, 600, r'$y = ' + '%.4f' % slope + 'x + ' + '%.4f' % intercept + '$', fontsize=16)  
    
18.     text(100, 520, r'$R^2 = '+ '%.4f' % (r * r) + '$', fontsize=16)
    
19.     text(700, 650, '(a)')
    
20. 
    
21. print 'Finished...'

fn = os.path.join('D:/KeyData/PMMUL/data/54500_PMMUL_DA.csv')<br /> if os.path.exists(fn):<br />     print fn<br />     tdata = readtable(fn, delimiter=',', format='%{yyyyMMdd}D%f%f%f%i%i%i')        <br />     pm10 = tdata['PM10']<br />     pm25 = tdata['PM2.5']<br />     slope, intercept, r, p, stderr, n = linregress(pm10, pm25)<br />     x = arange(0, 800, 100)<br />     y = polyval([slope, intercept], x)<br /> <br />     scatter(pm10, pm25, s=2, color='k', label='PM')<br />     plot(x, y, 'r-', linewidth=2)<br />     xlim(0, 800)<br />     ylim(0, 800)<br />     xlabel(r'$\rm{PM_{10}} \ (\mu g \ m^{-3})$')<br />     ylabel(r'$\rm{PM_{2.5}} \ (\mu g \ m^{-3})$') <br />     text(100, 600, r'$y = ' + '%.4f' % slope + 'x + ' + '%.4f' % intercept + '$', fontsize=16)  <br />     text(100, 520, r'$R^2 = '+ '%.4f' % (r * r) + '$', fontsize=16)<br />     text(700, 650, '(a)')<br /> <br /> print 'Finished...'

zihan6688

很不错的程序，谢谢您在MeteoInfo上的辛勤付出

bjydly

现下载备用着

zsmnantong

可能会用到，谢谢！

ldepn

运行的时候出现一个错误
错误出现在：tdata = readtable(fn, delimiter=',', format='%{yyyyMMddhh}D%f%f%f%i%i%i%i')      
提示如下：

E:/Baidudownload/MeteoInfo/script/beijingPM25.csv
Traceback (most recent call last):
  File "<iostream>", line 4, in <module>
  File "E:\Baidudownload\MeteoInfo\MeteoInfo(Java)\pylib\mipylib\midata.py", line 306, in readtable
    tdata = TableUtil.readASCIIFile(filename, delimiter, headerlines, format, encoding)
        at org.meteoinfo.data.TableUtil.readASCIIFile(TableUtil.java:78)
        at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
        at sun.reflect.NativeMethodAccessorImpl.invoke(Unknown Source)
        at sun.reflect.DelegatingMethodAccessorImpl.invoke(Unknown Source)
        at java.lang.reflect.Method.invoke(Unknown Source)

MeteoInfo

ldepn 发表于 2015-7-17 11:28
运行的时候出现一个错误
错误出现在：tdata = readtable(fn, delimiter=',', format='%{yyyyMMddhh}D%f%f% ...

你的数据可能和我用的示例数据格式不同。

readtable函数是将文本文件中的数据读入到一个二维表格中去，类似将数据导入到Excel中的一个sheet中。要求数据文件第一行是数据每列的列名（字段名），后面是一行行的数据，每行数据的列数应该一样。readtable函数中的delimiter参数指定了数据之间的分隔符，如果是空格分隔则不用指定此参数（是缺省值）。format参数指定了每列数据的数据类型，%f是浮点型（float），%s是字符型（string），%d是双精度浮点型（double），%i是整形（int），日期型比较特殊，需要用大括号，比如%{yyyyMMdd}D指定了日期是4位年2位月2位日组成的字符串。1楼用的示例数据格式如下：
Time,PM10,PM2.5,PM1,N_PM10,N_PM2.5,N_PM1
20060306,163.4,76.7,60.9,8,9,9
20060307,243.9,135.3,114.5,19,19,19
20060308,168.6,95.9,79.9,15,15,15
20060309,233.3,123.2,100.8,22,22,22
20060310,521.0,155.2,93.5,13,13,13
20060311,123.5,35.3,22.9,19,19,19
20060312,56.1,14.6,10.1,15,15,15
20060313,125.4,34.1,21.0,13,13,13


你需要根据自己的数据格式做相应的调整。读入数据后可以根据字段名取某一列数据。

ldepn

MeteoInfo 发表于 2015-7-17 12:31
你的数据可能和我用的示例数据格式不同。

readtable函数是将文本文件中的数据读入到一个二维表格中去 ...

如何读取数据的是时间有小时，“%{yyyyMMdd}”这个该怎么改？这个是什么语言的函数？在网上没有找到这个函数的使用格式

MeteoInfo

ldepn 发表于 2015-7-17 14:15
如何读取数据的是时间有小时，“%{yyyyMMdd}”这个该怎么改？这个是什么语言的函数？在网上没有找到这个 ...

你把你的数据前几行贴出来（类似楼上的帖子）看看。

ldepn

MeteoInfo 发表于 2015-7-17 14:31
你把你的数据前几行贴出来（类似楼上的帖子）看看。

Time        PM2.5        O3        NO        NO2        SO2        CO
2014/11/01 00:00        122.6        9.3375        51.57        40.5775        5.32        3.3215
2014/11/01 00:10        120.2        9.25        45.835        37.825        4.915        3.281
2014/11/01 00:20        105.98        9.31        43.095        38.62        4.805        3.263
2014/11/01 00:30        112.92        9.225        46.67        39.475        5.205        3.2605
2014/11/01 00:40        129        9.625        43.775        40.145        5.02        3.264
2014/11/01 00:50        122.24        9.505        44.405        40.92        4.99        3.2785
2014/11/01 01:00        120.96        9.285        43.7        39.28        4.855        3.2625
2014/11/01 01:10        105.6        9.08        44.085        36.675        4.555        3.2215
2014/11/01 01:20        117        9.46        44.215        39.64        4.645        3.268
2014/11/01 01:30        109.68        9.435        41.875        37.465        4.3475        3.234
2014/11/01 01:40        113.44        9.65        43.18        38.195        4.525        3.2145
2014/11/01 01:50        118.6        9.6        43.275        37.63        4.525        3.196
2014/11/01 02:00        103.32        9.3675        44.9175        36.555        4.4775        3.1885
2014/11/01 02:10        104.84        9.655        41.53        36.76        4.435        3.186
2014/11/01 02:20        112.04        9.4475        39.6425        35.0675        4.3425        3.1725
2014/11/01 02:30        106.48        9.785        39.9        35.525        4.31        3.1685
2014/11/01 02:40        104.58        10.4175        38.41        35.1475        4.365        3.087
2014/11/01 02:50        100.16        9.73        38.465        34.175        4.175        3.09
2014/11/01 03:00        113        10.805        36.755        33.52        4.32        2.9945
2014/11/01 03:10        97.94        10.385        35.6975        32.1275        4.195        2.971

这是在excel里面的格式

MeteoInfo

ldepn 发表于 2015-7-17 14:40
Time        PM2.5        O3        NO        NO2        SO2        CO
2014/11/01 00:00        122.6        9.3375        51.57        40.5775        5.32        3.3215
2014/11/01 0 ...

你可以另存为csv（逗号分隔）文件，然后用下面语句读入数据：

tdata = readtable(fn, delimiter=',', format='%{yyyy/MM/dd hh:mm}D%f%f%f%f%f%f')