notebook.community

Edit and run

如何使用Python讀寫檔案資料？

參考資料

準備工作：所需檔案下載

注意：python2 與 python3 下的 urllib有點不同, 下面修改成 python2 的環境, python3的語法會 comment 起來





In [1]:



    


# import urllib.request
import urllib

下面兩個 cell 分別演示如何使用 urllib來透過 url 抓取資料

  • SDSS的 fits binary table, 裡面存的是光譜
  • NOAA的太陽黑子預測, 純文字檔案

執行後會多出兩個檔案 "spec.fits"和 "astro_table.txt", 可以用以前熟悉的工具先看看資料的內容與格式, 這樣比較有 fu...





In [2]:



    


fitsurl="http://dr10.sdss3.org/sas/dr10/sdss/spectro/redux/26/spectra/0651/spec-0651-52141-0569.fits"
#urllib.request.urlretrieve(fitsurl,"spec.fits")
urllib.urlretrieve(fitsurl,"spec.fits")



    


















    
Out[2]:








('spec.fits', <httplib.HTTPMessage instance at 0x103fac440>)

















In [3]:



    


texturl="http://services.swpc.noaa.gov/text/predicted-sunspot-radio-flux.txt"
#urllib.request.urlretrieve(texturl, "astro_table.txt")
urllib.urlretrieve(texturl, "astro_table.txt")



    


















    
Out[3]:








('astro_table.txt', <httplib.HTTPMessage instance at 0x103fc0d88>)

Lab：請修改 url 為你自己感興趣的資料來源, 下載資料並存在自己命名的檔案中

範例：用Astropy 讀 FITS檔





In [14]:



    


from astropy.io import fits
hdulist = fits.open('spec.fits')
hdulist.info()



    


















    






Filename: spec.fits
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU     141   ()              
1    COADD       BinTableHDU     26   3829R x 8C   [E, E, E, J, J, E, E, E]   
2    SPECOBJ     BinTableHDU    262   1R x 126C    [6A, 4A, 16A, 23A, 16A, 8A, E, E, E, J, E, E, J, B, B, B, B, B, B, J, 22A, 19A, 19A, 22A, 19A, I, 3A, 3A, 1A, J, D, D, D, E, E, 19A, 8A, J, J, J, J, K, K, J, J, J, J, J, J, K, K, K, K, I, J, J, J, J, 5J, D, D, 6A, 21A, E, E, E, J, E, 24A, 10J, J, 10E, E, E, E, E, E, E, J, E, E, E, J, E, 5E, E, 10E, 10E, 10E, 5E, 5E, 5E, 5E, 5E, J, J, E, E, E, E, E, E, 25A, 21A, 10A, E, E, E, E, E, E, E, E, J, E, E, J, 1A, 1A, E, E, J, J, 1A, 5E, 5E]   
3    SPZLINE     BinTableHDU     48   29R x 19C    [J, J, J, 13A, D, E, E, E, E, E, E, E, E, E, E, J, J, E, E]   
4    B2-00010589-00010594-00010595  BinTableHDU    146   2047R x 7C   [E, E, E, J, E, E, E]   
5    B2-00010590-00010594-00010595  BinTableHDU    146   2047R x 7C   [E, E, E, J, E, E, E]   
6    B2-00010591-00010594-00010595  BinTableHDU    146   2047R x 7C   [E, E, E, J, E, E, E]   
7    B2-00010592-00010594-00010595  BinTableHDU    146   2047R x 7C   [E, E, E, J, E, E, E]   
8    R2-00010589-00010594-00010595  BinTableHDU    146   2044R x 7C   [E, E, E, J, E, E, E]   
9    R2-00010590-00010594-00010595  BinTableHDU    146   2044R x 7C   [E, E, E, J, E, E, E]   
10   R2-00010591-00010594-00010595  BinTableHDU    146   2044R x 7C   [E, E, E, J, E, E, E]   
11   R2-00010592-00010594-00010595  BinTableHDU    146   2044R x 7C   [E, E, E, J, E, E, E]

hdulist 顧名思義, 就是 a list of Header/Data Unit (廢話), 關於 FITS檔案的格式, 請參閱 http://fits.gsfc.nasa.gov/fits_primer.html 裡面有詳盡的說明

hdulist.info() 會把每個 element的摘要列出來, 以上面的例子來說, 就是讀進來的 spec.fits 他有 12個 units, 第 0個就是個純的 header, 第 1,2,3個 unit應該是一些 meta data, 真正的科學資料是從第 4個 unit開始, 格式為 binary table, 接下來我們就用第四個 unit 來作演示

通常在不太了解新的物件時, dir() 是我們的好朋友, 方便我們猜測可以對這個陌生的物件做些什麼事情。 所以, 就 dir()吧





In [9]:



    


dir(hdulist[4])



    


















    
Out[9]:








['_EXCLUDE',
 '_MASK',
 '__class__',
 '__delattr__',
 '__dict__',
 '__doc__',
 '__format__',
 '__getattribute__',
 '__hash__',
 '__init__',
 '__module__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 '_buffer',
 '_calculate_checksum',
 '_calculate_datasum',
 '_calculate_datasum_with_heap',
 '_char_encode',
 '_clear_table_keywords',
 '_close',
 '_columns_type',
 '_compute_checksum',
 '_compute_hdu_checksum',
 '_datLoc',
 '_datSpan',
 '_data_loaded',
 '_data_needs_rescale',
 '_data_offset',
 '_data_replaced',
 '_data_size',
 '_data_type',
 '_default_name',
 '_dump_coldefs',
 '_dump_data',
 '_encode_byte',
 '_ext_comment',
 '_extension',
 '_file',
 '_get_raw_data',
 '_get_tbdata',
 '_get_timestamp',
 '_has_data',
 '_hdrLoc',
 '_hdu_registry',
 '_header',
 '_header_offset',
 '_init_tbdata',
 '_load_coldefs',
 '_load_data',
 '_manages_own_heap',
 '_new',
 '_nrows',
 '_output_checksum',
 '_padding_byte',
 '_populate_table_keywords',
 '_postwriteto',
 '_prewriteto',
 '_readfrom_internal',
 '_standard',
 '_summary',
 '_tdump_file_format',
 '_theap',
 '_uint',
 '_update_checksum',
 '_update_column_added',
 '_update_column_attribute_changed',
 '_update_column_removed',
 '_update_uint_scale_keywords',
 '_verify',
 '_verify_checksum_datasum',
 '_writedata',
 '_writedata_by_row',
 '_writedata_direct_copy',
 '_writedata_internal',
 '_writeheader',
 '_writeto',
 'add_checksum',
 'add_datasum',
 'columns',
 'copy',
 'data',
 'dump',
 'filebytes',
 'fileinfo',
 'from_columns',
 'fromstring',
 'header',
 'is_image',
 'level',
 'load',
 'match_header',
 'name',
 'readfrom',
 'register_hdu',
 'req_cards',
 'run_option',
 'size',
 'tcreate',
 'tdump',
 'unregister_hdu',
 'update',
 'update_ext_name',
 'update_ext_version',
 'ver',
 'verify',
 'verify_checksum',
 'verify_datasum',
 'writeto']

雖然落落長一串, 希望你的眼睛沒有什麼業障, 可以發現有一些有趣的屬性/方法, 像是 data, dump, writeto... 等等 接下來就繼續冒險前進, 用點點 (.) 來探索新世界吧

所以就 .data





In [15]:



    


hdulist[4].data



    


















    
Out[15]:








FITS_rec([(17.107233, 3.790184, 0.00034346466, 0, 0.77565843, 16.30872, 0.42700836),
       (23.673752, 3.7900956, 0.00058561546, 0, 0.77561331, 17.235472, 0.41987428),
       (28.730856, 3.7900071, 0.00084290426, 0, 0.77556837, 17.819456, 0.41240782),
       ...,
       (31.470488, 3.5814037, 0.012394244, 0, 1.0769527, 12.84608, 0.59983635),
       (30.73415, 3.5812922, 0.011396203, 0, 1.0772977, 12.612031, 0.60809124),
       (39.182308, 3.5811803, 0.0095504085, 0, 1.0776428, 12.544607, 0.62002754)], 
      dtype=(numpy.record, [('flux', '>f4'), ('loglam', '>f4'), ('ivar', '>f4'), ('mask', '>i4'), ('wdisp', '>f4'), ('sky', '>f4'), ('calib', '>f4')]))

你看看你, data這不就來了嗎~~~

好, 看來這像是個 numpy 的 record array (or structured array), 格式長這樣 : ([(v1, v2, v3...), (), .... ()], dtype=(type1, type2, type3...)), 不難猜測 v1 對應 type1... 以此類推

  • 所以, 17.107是 flux, 3.79是 loglam.... 等等等等

我想 loglam是 lambda 的 log這種事應該猜的到齁？還猜不到的話, 恭喜你, 還中毒未深~~~ 要抓出 flux 或任意 type的方法也很直覺





In [16]:



    


hdulist[4].data["flux"]



    


















    
Out[16]:








array([ 17.10723305,  23.67375183,  28.73085594, ...,  31.47048759,
        30.73414993,  39.1823082 ], dtype=float32)

好滴... 這個演示就先告個段落, 記得要把打開的東西 (hdulist) 關起來, 養成好習慣以後才不會很難改





In [11]:



    


hdulist.close()

範例：用Numpy讀寫檔案





In [6]:



    


# 尚需要更合適的檔案來當範例



    











In [21]:



    


import numpy as np



    











In [ ]:



    


讀進的array = np.loadtxt("檔名")



    











In [ ]:



    


# 對讀進的array做一些運算後存成新array，然後將新array寫到一個新檔案中 (代補)



    











In [ ]:



    


np.savetxt('新檔名', 新array)

範例：用Astropy讀寫 ASCII 檔案

# 尚需要更合適的檔案來當範例
# 更新使用 NOAA的太陽黑子預報 http://services.swpc.noaa.gov/text/predicted-sunspot-radio-flux.txt

astropy 的 ascii 提供相當完整的 ascii 讀取方案, 除了基本的格式外, 還支援 csv, cds (vizier), latex, DAOPHOT, Sextractor... 天文領域中常見的 ascii 表格大概都可以無痛讀取, 不需要自己花時間解析

用法也很直覺, 如下:





In [18]:



    


from astropy.io import ascii
#data = ascii.read("astro_table.txt")
data = ascii.read("astro_table.txt",data_start=2)  
print(data)



    


















    






col1 col2 col3 col4 col5  col6  col7 col8
---- ---- ---- ---- ---- ----- ----- ----
2016    1 33.7 34.7 32.7 100.3 101.3 99.3
2016    2 33.9 35.9 31.9  99.5 100.5 98.5
2016    3 33.6 36.6 30.6  98.7 100.7 96.7
2016    4 33.2 38.2 28.2  97.9 100.9 94.9
2016    5 33.1 38.1 28.1  96.7 100.7 92.7
2016    6 33.0 39.0 27.0  95.0  99.0 91.0
2016    7 33.0 40.0 26.0  93.5  98.5 88.5
2016    8 32.8 39.8 25.8  92.2  98.2 86.2
2016    9 32.6 40.6 24.6  91.3  98.3 84.3
2016   10 32.7 41.7 23.7  90.7  98.7 82.7
 ...  ...  ...  ...  ...   ...   ...  ...
2019    2  8.3 18.3  0.0  66.3  75.3 60.0
2019    3  7.8 17.8  0.0  65.9  74.9 60.0
2019    4  7.3 17.3  0.0  65.4  74.4 60.0
2019    5  6.8 16.8  0.0  65.0  74.0 60.0
2019    6  6.4 16.4  0.0  64.5  73.5 60.0
2019    7  5.9 15.9  0.0  64.1  73.1 60.0
2019    8  5.5 15.5  0.0  63.8  72.8 60.0
2019    9  5.1 15.1  0.0  63.4  72.4 60.0
2019   10  4.8 14.8  0.0  63.1  72.1 60.0
2019   11  4.4 14.4  0.0  62.8  71.8 60.0
2019   12  4.1 14.1  0.0  62.5  71.5 60.0
Length = 48 rows

read() 出來的是個 table 物件, dir()的樂趣我就留給各位, 下面直接舉兩個常用的例子:

  • 抽取特定的 column(s)
  • 使用 where() 來抓取符合條件的 rows




In [19]:



    


print(data["col5","col1"])



    


















    






col5 col1
---- ----
32.7 2016
31.9 2016
30.6 2016
28.2 2016
28.1 2016
27.0 2016
26.0 2016
25.8 2016
24.6 2016
23.7 2016
 ...  ...
 0.0 2019
 0.0 2019
 0.0 2019
 0.0 2019
 0.0 2019
 0.0 2019
 0.0 2019
 0.0 2019
 0.0 2019
 0.0 2019
 0.0 2019
Length = 48 rows

















In [23]:



    


d = np.where((data["col1"] >= 2017) & (data["col1"] <= 2019) & (data["col2"] >=6) & (data["col2"] <=9))
print(d)
print(data[d]["col1","col2","col3"])



    


















    






(array([17, 18, 19, 20, 29, 30, 31, 32, 41, 42, 43, 44]),)
col1 col2 col3
---- ---- ----
2017    6 26.1
2017    7 24.9
2017    8 23.7
2017    9 22.5
2018    6 13.7
2018    7 12.9
2018    8 12.2
2018    9 11.5
2019    6  6.4
2019    7  5.9
2019    8  5.5
2019    9  5.1

















In [ ]:



    


dir(data)



    











In [ ]:



    


data.write("out", format="ascii")

Content source: Astrohackers-TW/IANCUPythonAdventure

Similar notebooks:

notebook.community | gallery | about