iT邦幫忙

2022 iThome 鐵人賽

DAY 7
0
AI & Data

氣象食材系列 第 7

[ Day 7] 基本觀測資料-雷達資料 (xml)

  • 分享至 

  • xImage
  •  

今天介紹雷達資料,其實精準地來說,是要看該雷達的功用,才能知道可以得到的資料。
譬如氣象雷達的功用,是透過發射電磁波,當電磁波遇到空氣中有足夠大的雨滴或冰晶時,就能夠測得相關資訊。

而最常用的相關資訊,就是雷達回波。故先介紹讀取雷達回波。

先至氣象局開放資料下載雷達回波資料,如下圖

https://ithelp.ithome.com.tw/upload/images/20220921/201509237c063omzMn.jpg

因為之前已經介紹過如何讀取xml,故下載xml檔案格式,下載之後應該會得到O-A0059-001.xml檔案

使用網頁瀏覽器開啟檔案,如下圖
https://ithelp.ithome.com.tw/upload/images/20220921/20150923VHmZ6I3Aoo.jpg

可以觀察到所需的標記不多,所以就直接iter想要的標記。
處理的方式如下

import xml.etree.ElementTree as ET
import numpy as np
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

trees = ET.parse('O-A0059-001.xml')
roots = trees.getroot()

pN = "{urn:cwb:gov:tw:cwbcommon:0.1}parameterName"
rN = "{urn:cwb:gov:tw:cwbcommon:0.1}radarName"
pV = "{urn:cwb:gov:tw:cwbcommon:0.1}parameterValue"

parval = []
for hdvalue in roots.iter(pV):
    parval.append(hdvalue.text)

dictparm = dict()
cnt = 0
for hd in roots.iter(pN):
    if hd.text == '整合雷達名':
        for hdvalue in roots.iter(rN):
            dictparm.update({hd.text:hdvalue.text})
    else:
        dictparm.update({hd.text:parval[cnt]})
        cnt +=1
print(dictparm)

ac = []
for qq in roots.iter("{urn:cwb:gov:tw:cwbcommon:0.1}content"):
    ac = qq.text
    
reflect = ac.split(",")
reflect = np.array([float(x) for x in reflect]).reshape(881,921)

fig = plt.figure(figsize=(12,8))
axs = plt.axes(projection=crs.PlateCarree())

tw_shp = ShapelyFeature(shpreader.Reader("C:\workplace\ithome\/tw_shp\COUNTY_MOI_1090820.shp").geometries(), crs.PlateCarree(),facecolor="none",edgecolor='k')
lon,lat = np.meshgrid(np.linspace(115.0, 115.0+0.0125*920,921), np.linspace(18.0, 18.0+0.0125*880, 881))
cwb_data = ["#F2F2F2","white","#07FDFD","#0695FD", "#0203F9","#00FF00","#00C800","#019500","#FEFD02","#FEC801","#FD7A00","#FB0100","#C70100","#950100","#FA03FA","#9800F6"]
cmaps = mcolors.ListedColormap(cwb_data,'radar')
colorlevel = [-1000, -900, 0 ,5,10,15,20,25,30,35,40,45,50,55,60,65,70]
norms = mcolors.BoundaryNorm(colorlevel, cmaps.N)
axs.gridlines(draw_labels=True)
pcm = axs.pcolormesh(lon,lat,reflect,cmap = cmaps, norm=norms,transform=crs.PlateCarree())
axs.add_feature(tw_shp)
axs.set_extent([lon[0,0],lon[-1,-1],lat[0,0],lat[-1,-1]])
plt.title(dictparm["時間"])
#故意做一個colorbar 
cwb_data = ["#07FDFD","#0695FD", "#0203F9","#00FF00","#00C800","#019500","#FEFD02","#FEC801","#FD7A00","#FB0100","#C70100","#950100","#FA03FA","#9800F6"]
colorlevel = [ 0 ,5,10,15,20,25,30,35,40,45,50,55,60,65,70]
cmaps = mcolors.ListedColormap(cwb_data,'radar')
norms = mcolors.BoundaryNorm(colorlevel, cmaps.N)
plt.colorbar(mpl.cm.ScalarMappable(norm=norms, cmap=cmaps),pad=0.1)

這樣就可以讀取雷達回波跟視覺化囉

視覺化的圖如下
https://ithelp.ithome.com.tw/upload/images/20220921/201509235fk9Cdcnbk.png

驗證一下,取氣象局相同時間的雷達回波圖
https://ithelp.ithome.com.tw/upload/images/20220921/20150923UMwquI3MXk.png

看起來應該是對的喔~


上一篇
[ Day 6] 基本觀測資料-氣象站資料4 (xml)
下一篇
[Day 8] 基本觀測資料-過去1小時雷達定量降雨估計 (xml)
系列文
氣象食材30
圖片
  直播研討會
圖片
{{ item.channelVendor }} {{ item.webinarstarted }} |
{{ formatDate(item.duration) }}
直播中

尚未有邦友留言

立即登入留言