---
# System prepended metadata

title: Create StationXML with ObsPy
tags: [Python]

---

---
GA: G-44K4TH870N
---
> The purpose of the FDSN StationXML schema is to define an XML representation of the most important and commonly used metadata structures to describe digital seismological data, in particular for the purposes of data use.
### StationXML
隨著地震學研究逐漸和資訊工程接軌，許多傳統的資料格式也逐漸統一朝向更泛用和容易擴充的格式。其中FDSN制定的[StationXML](https://docs.fdsn.org/projects/stationxml/en/latest/index.html)格式就被用來儲存地震測站的各種資訊。

在StationXML之前，儲存地震資料通用的格式是[SEED](https://ds.iris.edu/ds/nodes/dmc/data/formats/seed/) （Standard for the Exchange of Earthquake Data），而沒有包含地震資料而只有測站資訊的格式則稱為Dataless SEED（或簡稱 Dataless）。因此StationXML的設計主要還是繼承了SEED原有的欄位，以[XML](https://zh.wikipedia.org/zh-tw/XML)（一種階層清單格式，因為通用性很高所以被廣泛使用，最常見的即為網頁所使用的HTML）的格式儲存，更容易被現在的工具讀寫。

StationXML被稱為是一個清單（Inventory），其實可以包含多個測網（Network）或多個測站。

![image](https://hackmd.io/_uploads/SJljmAdwR.png)

可以看到，一個Inventory下可以有複數個Network，每個Network可以有多個Station，每個Station又可以有多個Channel（例如三軸地震儀的ZNE三頻道）。而每個階層的詳細資料就並列儲存於該階層中。

BATS地震站TWGB的Inventory範例
```xml!
<?xml version='1.0' encoding='UTF-8'?>
<FDSNStationXML xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.2">
  <Source>SeisComP</Source>
  <Sender>CWB</Sender>
  <Module/>
  <ModuleURI/>
  <Created>2023-08-09T02:17:01.898734Z</Created>
  <Network code="TW" startDate="1980-01-01T00:00:00.000000Z">
    <Description>Seismic Array in Taiwan</Description>
    <Station code="TWGB" startDate="1996-03-08T00:00:00.000000Z">
      <Latitude unit="DEGREES">22.8177</Latitude>
      <Longitude unit="DEGREES">121.0799</Longitude>
      <Elevation>289.9</Elevation>
      <Site>
        <Name>Taitung</Name>
      </Site>
      <CreationDate>1996-03-08T00:00:00.000000Z</CreationDate>
      <Channel code="BHE" startDate="1996-03-08T00:00:00.000000Z" endDate="1997-07-26T02:29:00.000000Z" locationCode="">
```
### ObsPy的Inventory

我們使用Python下最大的地震學套件ObsPy來製作一份StationXML：
https://docs.obspy.org/tutorial/code_snippets/stationxml_file_from_scratch.html

在Obspy.core.inventory中有對應到StationXML不同階層的Class，也就是Networks-Stations-Channels這幾個型別。所以實作上就是個別建立這些型別的變數，並把資料儲存在裡面，最後再一層一層Append到最上層。

### 把儀器響應加進Inventory

再來就要根據手上現有的資料格式再設計程式去讀取資料並寫進變數。例如我有純粹的PoleZero儀器響應資料，但是沒有測站的經緯度，那我就需要分別讀取資料。

一個典型的PoleZero響應要包含至少三個參數，Poles、Zeros、Constant(Gain)。通常看起來會像這樣：

```
ZEROS	4
	-4.341000e+02	+0.000000e+00
	+0.000000e+00	+0.000000e+00
	+0.000000e+00	+0.000000e+00
	+0.000000e+00	+0.000000e+00
POLES	7
	-3.691000e-02	+3.712000e-02
	-3.691000e-02	-3.712000e-02
	-3.712000e+02	+0.000000e+00
	-3.739000e+02	+4.755000e+02
	-3.739000e+02	-4.755000e+02
	-5.884000e+02	+1.508000e+03
	-5.884000e+02	-1.508000e+03
CONSTANT	2.5713719863919994e+20	
```

那我們就可以寫一個簡單的函數來讀它：

```python=
def read_pz_file(filename):
    zeros = []
    poles = []
    constant = None
    with open(filename, "r") as file:
        lines = file.readlines()
        i = 0
        while i < len(lines):
            line = lines[i].strip()
            if line.startswith("ZEROS"):
                num_zeros = int(line.split()[1])
                for j in range(num_zeros):
                    i += 1
                    zero_parts = lines[i].split()
                    zeros.append(complex(float(zero_parts[0]), 
                                         float(zero_parts[1])))
            elif line.startswith("POLES"):
                num_poles = int(line.split()[1])
                for j in range(num_poles):
                    i += 1
                    pole_parts = lines[i].split()
                    poles.append(complex(float(pole_parts[0]),
                                         float(pole_parts[1])))
            elif line.startswith("CONSTANT"):
                constant = float(line.split()[1])
            i += 1
    return zeros, poles, constant
```
接著再用ObsPy的這個Class method，Response.from_paz()把這三個參數變成Response格式，加到這個channel中，再把channel給加進測站的channel清單。
```python=
this_channel_response = Response.from_paz(zeros, poles, constant)
channel.response = this_channel_response
station.channels.append(channel)
```
如此一來我們就把儀器響應的資料從polezero檔匯入測站的Inventory囉！