earth套件貼上CAMS臭氧濃度

Table of contents

背景
- CAMS與GFS數據之間需合併的項目
- 色階的選擇
數據下載與轉換
earth系統新增臭氧之讀取繪圖功能
- html程式
- products.js程式
結果討論

背景

雖然臭氧是光化煙霧的重要指標，也是本土石化業與汽機車污染所影響的空品項目，然而earth.nullschool或是windy目前都沒有貼上臭氧濃度之實例，原因不明，可能並不是很多人這麼瞭解臭氧的指標意義吧。
ECWMF的CAMS全球空品預報有這個項目，解析度為0.4度每12小時更新，貼在背景1度的GFS風場，檔案小、反應快、有其便利性。
- CAMS有每半年更新的再分析數據、有近實時再分析、也有空品預報數據，這三者雖然精確性以最後者最低，但時效性卻最高。
- GFS雖然也有提供臭氧的混合比(var_O3MR)，但只有1000mb等壓面值，並沒有地面值，因此也無法應用。
- CAMS與GFS數據之間需合併的項目如表所示

CAMS與GFS數據之間需合併的項目

項目	CAMS	GFS	說明
空間解析度	0.4度	1度	2者解析度不能整除，還是需要進行內插
經度0度	沒有值	有值	前者必須另行自360度取值
更新頻率	12小時	6小時	需由不同的crontab腳本執行下載與處理
預報時間(leadtime_hour)	3小時	3小時	相同
官方圖面	哥白尼網站	NWS	後者沒有展示美國以外地區預報結果

把下載與應用的流程自動化後，不僅可以提供檢視目前的氣流與光化煙霧跨境傳輸現象，也可以使高解析度空品即時模擬向前推動一大步。

色階的選擇

earth.nullschoolchem的色階多應用在原生性污染物，其特性為少數高值，大多位置濃度均不高。用在濃度差異不大的空品項目圖面顏色將會太少。
參考earth.nullschool的Particulates展現方式，基本上為一彩虹色階。
雖然背景臭氧濃度高低差異不是很大，但為了將來模式模擬解析度如果提高，會出現局部超高/超低濃度，彩虹色階會有較佳的表現。
也能接近環保署官方網站的展示方式。

數據下載與轉換

數據量雖然不大，但還需在CAMS排隊等候下載，會需要一些時間(~20分鐘)。
數據下載
詳歐洲中期天氣預報中心再分析數據之下載
- 對話網頁：https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-atmospheric-composition-forecasts?tab=form
以下為2022-08-01下載程式範例

#kuang@master /nas1/ecmwf/CAMS/CAMS_global_atmospheric_composition_forecasts/2022
#$ cat get_O3.py
import cdsapi

c = cdsapi.Client()

c.retrieve(
    'cams-global-atmospheric-composition-forecasts',
    {
        'date': '2022-08-01/2022-08-01',
        'type': 'forecast',
        'format': 'grib',
        'time': '12:00',
        'pressure_level': '1000',
        'variable': [ 'ozone', ],
        'leadtime_hour': '0',
    },
    'ozone_globe.grib')

有關level

CAMS高度定義有2種

Pressure levels(mb)
- 1000/950/925/900/850/800/700/600/500/400/300/250/200/150/100/70/50/30/20/10/7/5/3/2/1
- 共25層
- 由於定壓層高度會隨時間地點改變，並不用在空品模擬。無法使用。
Model levels
- 1(top)to 137(surface),
- which are described at https://confluence.ecmwf.int/display/UDOC/L137+model+level+definitions.
- Before 9 July 2019, the vertical levels were 60, which are described at https://confluence.ecmwf.int/display/UDOC/L60+model+level+definitions.
從官網取得table.csv後進行處理如下
- bisect即使是bisect_left還是取不到地面，故將其減1。

df=read_csv('/nas1/ecmwf/CAMS/CAMS_global_atmospheric_composition_forecasts/2022/heights.csv')
h='Geometric Altitude [m]'
df=df.sort_values(h)
h137=list(df[h])
n137=list(df['n'])
kk=[bisect(h137,zz[k])-1 for k in range(24)]
[str(n137[kk[k]]) for k in range(24)]

代入get.py中 kk= ['137', '135', '133', '129', '125', '122', '120', '117', '114', '112', '110', '107', '105', '101', '96', '92', '87', '83', '78', '73', '67', '61', '56', '51']

數據轉換

grb2json.py@github
CAMS檔案內容之讀取
- CAMSgrib2檔案之y軸順序：自北向南(與GFS相同)
- 因需與座標系統搭配以進行內插，此處將其翻轉以避免錯誤。

fname=sys.argv[1]
grbs = pygrib.open(fname)
uv=['ozone', ]
atbs={'ozone': 'GEMS Ozone',}
for a in set(atbs):
  grb = grbs.select(name=atbs[a])
  cmd=a+'=grb[0].values'
  exec(cmd)
ozone=np.flip(ozone,axis=0)

使用griddata進行內插
- CAMS網格系統：xyc
- GFS網格系統：x1,y1
- GFSjson檔案之y軸順序：自北向南(與CAMS相同)

for ir in range(nr):
  c = np.array([ozone[idx[0][i], idx[1][i]] for i in range(mp)])*2.e9
  zz = griddata(xyc, c[:], (x1, y1), method='cubic')
  gfs[ir]['data']=list(np.flip(np.where(zz!=zz,0,zz),axis=0).flatten())

執行結果放在/Users/Data/javascripts/D3js/earth/public/data/weather/current以利js程式讀取
結果檔名：current-ozone-surface-level-gfs-1.0.json

earth系統新增臭氧之讀取繪圖功能

html程式

在下拉對話框增加貼上O3數據
- 直接在氣象數據後面接上空品項目
- 並不像nullschool將污染項目分類，以減少程式修改。

#kuang@114-32-164-198 /Users/Data/javascripts/D3js/earth
#$ diff /Users/Data/javascripts/D3js/earthGFS/public/index.html public/index.html
104c104,105
<                 class="text-button" id="overlay-mean_sea_level_pressure" title="Mean Sea Level Pressure">MSLP</span>
---
>                 class="text-button" id="overlay-mean_sea_level_pressure" title="Mean Sea Level Pressure">MSLP</span> – <span
>                 class="text-button" id="overlay-ozone" title="GEMS Ozone">O3</span>

products.js程式

程式參考溫度(temp)段落
色階也參考溫度的設定
- 取消粉紅及青色(cyan)，顏色太過顯眼、彩度不足、也不是一般彩虹色階項目

#$ diff /Users/Data/javascripts/D3js/earthGFS/public/libs/earth/1.0.0/products.js public/libs/earth/1.0.0/products.js
401a408,453
>
>         "ozone": {
>             matches: _.matches({param: "wind", overlayType: "ozone"}),
>             create: function(attr) {
>                 return buildProduct({
>                     field: "scalar",
>                     type: "ozone",
>                     description: localize({
>                         name: {en: "GEMS Ozone", ja: "臭氧"},
>                         qualifier: {en: " @ " + describeSurface(attr), ja: " @ " + describeSurfaceJa(attr)}
>                     }),
>                     paths: [gfs1p0degPath(attr, "ozone", attr.surface, attr.level)],
>                     date: gfsDate(attr),
>                     builder: function(file) {
>                         var record = file[0], data = record.data;
>                         return {
>                             header: record.header,
>                             interpolate: bilinearInterpolateScalar,
>                             data: function(i) {
>                                 return data[i];
>                             }
>                         }
>                     },
>                     units: [
>                         {label: "ppb", conversion: function(x) { return x; }, precision: 3}
>                     ],
>                     scale: {
>                         bounds: [0, 500],
>                         gradient:
>                             µ.segmentedColorScale([
>                             [  0,  [37, 4, 42]],
>                             [ 30,     [41, 10, 130]],//purple blue
>                             [ 60,  [21, 84, 187]],   //blue
>                             [ 90,  [24, 132, 14]],   //green
>                             [120,  [247, 251, 59]], //yellow
>                             [150,  [235, 167, 21]],//
>                             [180,  [230, 71, 39]],
>                             [300,  [88, 27, 67]],
>                             [500,  [81, 40, 40]],//brown
>                             ])
>                     }
>                 });
>             }
>         },
>                             //[ 90,  [192, 37, 149]],  //pink
>                             //[120, [70, 215, 215]],  //cyan

結果討論

實例網址：http://sinotec24.com:8080
數據提交earth系統之圖面


earth貼上臭氧濃度之色階應用(2022080112)

哥白尼官網之圖面(下拉選項，不能放大縮小)

| | |:–:| | 同一時間哥白尼官網之圖面|

前者高濃度較為明顯、後者限於13個濃度層級，數字對照較為清晰。
差異可能來源：雖同為地面濃度，但前者未校正空氣密度，在內陸高原地區造成高估