Day 7 介紹雷達迴波資料，而雷達回波愈強，代表觀測到的雨滴或冰晶愈大，可能就會有明顯的降水。
所以過去氣象學家發展了雷達回波與降水強度的關係式，稱為ZR關係式。
透過ZR關係式推求的降水，即稱為雷達定量降雨估計。

氣象局的開放資料有提供這項資料

一樣下載xml檔案，使用網頁瀏覽器開啟

可以看到需要的資料是在content標記。如果是用xml.etree.ElementTree，只要先將根物件取出，接者iter即可，詳情可以參考Day 7的方式。

這邊介紹xmltodict套件處理xml檔案，安裝方式如下，
conda install -c conda-forge xmltodict
或
pip3 install xmltodict

安裝完成後，就開始使用吧

先讀取xml檔案(import numpy 目的是將非陣列資料變為陣列物件)

import xmltodict as xdict
import numpy as np

with open("O-B0045-001.xml", "r", encoding="utf-8") as fn:
    qpedict = xdict.parse(fn.read())

讀取之後，會得到複合dict物件，所以就可以一層一層拆解。以下拆解到第4層的key值

for key in qpedict.keys():
    print(1,key)
    for subkey in qpedict[key].keys():
        print(2,subkey)
        if isinstance(qpedict[key][subkey], dict):
            for subsubkey in qpedict[key][subkey].keys():
                print(3,subsubkey)
                if isinstance(qpedict[key][subkey][subsubkey], dict):
                    for subsubsubkey in qpedict[key][subkey][subsubkey].keys():
                        print(4,subsubsubkey)

結果如下
1 cwbopendata
2 @xmlns
2 identifier
2 sender
2 sent
2 status
2 msgType
2 scope
2 dataid
2 source
2 dataset
3 datasetInfo
4 datasetDescription
4 parameterSet
3 contents
4 contentDescription
4 content

可以知道複合dict物件的第一個key為cwbopendata，先取cwbopendata後，第二層的key值就會有@xmlns、identifier、sender、sent、status、msgType、scope、dataid、source、dataset，然後呢，可以看到dataset層還有第3層，分別為datasetInfo跟contents，剩下的以此類推。那我們要怎麼取值呢？其實就跟用dict物件一樣，今天想要看第三層datasetInfo的內容，就要用qpedict["cwbopendata"]["dataset"]["datasetInfo"]，這樣就可以取得。在這樣的基礎下，開始取我們需要的資料及資訊囉。

llcorlon, llcorlat = qpedict["cwbopendata"]["dataset"]["datasetInfo"]["parameterSet"]['parameter'][0]["parameterValue"].split(",")
llcorlon = float(llcorlon)
llcorlat = float(llcorlat)
res = float(qpedict["cwbopendata"]["dataset"]["datasetInfo"]["parameterSet"]['parameter'][1]["parameterValue"])
obstime = qpedict["cwbopendata"]["dataset"]["datasetInfo"]["parameterSet"]['parameter'][2]["parameterValue"]
gridx,gridy = qpedict["cwbopendata"]["dataset"]["datasetInfo"]["parameterSet"]['parameter'][3]["parameterValue"].split("*")
gridx = int(gridx)
gridy = int(gridy)

#因為取出來的資料都是str，故有些需要轉換成int或float。

lon, lat = np.meshgrid(np.linspace(llcorlon, llcorlon+res*(gridx-1), gridx),
                       np.linspace(llcorlat, llcorlat+res*(gridy-1), gridy))

#取降雨
qpestr = qpedict["cwbopendata"]["dataset"]["contents"]["content"].split(",")
qpe =np.array([ float(x) for x in rainstr]).reshape(gridy,gridx)

視覺化如下

import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.colors as mcolors
import cartopy.crs as crs
import cartopy.io.shapereader as shpreader
from cartopy.feature import ShapelyFeature

colorlevel = [0,1,2,6,10,15,20,30,40,50,70,90,110,130,150,200,300,400]
cwb_data=['None','#9BFFFF','#00CFFF','#0198FF','#0165FF','#309901','#32FF00','#F8FF00','#FFCB00',\
               '#FF9A00','#FA0300','#CC0003', '#A00000','#98009A','#C304CC','#F805F3','#FECBFF']
cmaps = mcolors.ListedColormap(cwb_data,'precipitation')
norms = mcolors.BoundaryNorm(colorlevel, cmaps.N)
fig = plt.figure(figsize=(12,8))
axs = plt.axes(projection=crs.PlateCarree())
tw_shp = ShapelyFeature(shpreader.Reader("D:\ithome\/tw_shp\COUNTY_MOI_1090820.shp").geometries(), crs.PlateCarree(),facecolor="none",edgecolor='k')
axs.gridlines(draw_labels=True)
ctf = axs.contourf(lon,lat,qpe,cmap = cmaps, norm=norms,levels=colorlevel,transform=crs.PlateCarree())
axs.add_feature(tw_shp)
plt.title(obstime)
plt.colorbar(ctf,ticks=colorlevel[1:len(colorlevel)-1],pad=0.07)