{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8146c21d",
   "metadata": {},
   "source": [
    "## 背景\n",
    "- 高空點源的**時變係數**骨幹是CEMS數據，然而同一工廠無數據、鄰近工業區其他廠無數據者，亦會參考CEMS設定其**時變係數**。\n",
    "- 排放量整體處理原則參見[處理程序總綱](https://sinotec2.github.io/Focus-on-Air-Quality/EmisProc/#處理程序總綱)、針對[點源之處理](https://sinotec2.github.io/Focus-on-Air-Quality/EmisProc/ptse/)及[龐大`.dbf`檔案之讀取](https://sinotec2.github.io/Focus-on-Air-Quality/EmisProc/dbf2csv.py/)，為此處之前處理。程式也會呼叫到[ptse_sub](https://sinotec2.github.io/Focus-on-Air-Quality/EmisProc/ptse/ptse_sub/)中的副程式\n",
    "\n",
    "## 程式說明\n",
    "\n",
    "### 排放與CEMS資料檔之讀取及準備\n",
    "- 引用模組\n",
    "  - 程式用到[ptse_sub](https://sinotec2.github.io/Focus-on-Air-Quality/EmisProc/ptse/ptse_sub/)中的副程式`CORRECT`, `add_PMS`, `check_nan`, `check_landsea`, `FillNan`, `WGS_TWD`, `Elev_YPM`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "59b647e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas import *\n",
    "import numpy as np\n",
    "import os, sys, subprocess\n",
    "import netCDF4\n",
    "import twd97\n",
    "import datetime\n",
    "from calendar import monthrange\n",
    "from scipy.io import FortranFile\n",
    "\n",
    "from ptse_sub import CORRECT, add_PMS, check_nan, check_landsea, FillNan, WGS_TWD, Elev_YPM\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19dd1019",
   "metadata": {},
   "source": [
    "- 從工作目錄讀取teds版本與年代"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3b458df5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Main\n",
    "P=subprocess.check_output('pwd',shell=True).decode('utf8').strip('\\n')+'/'\n",
    "teds=int(P.split('/')[3][-2:])\n",
    "yr=2016+(teds-10)*3\n",
    "ndays=365\n",
    "if yr%4==0:ndays=366\n",
    "s365=set([i*24 for i in range(ndays)])\n",
    "nhrs=ndays*24"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cca6fd07",
   "metadata": {},
   "source": [
    "- 檔案讀取與品質確認\n",
    "  - 此個案將所有點源資料庫的數據都以「高空」方式處理(cutting height of stacks`Hs=0`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3f1a422f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "set()\n",
      "TSP_EMI 0\n",
      "PM_EMI 0\n",
      "PM25_EMI 0\n",
      "SOX_EMI 0\n",
      "NOX_EMI 0\n",
      "THC_EMI 0\n",
      "NMHC_EMI 0\n",
      "CO_EMI 0\n",
      "PB_EMI 0\n"
     ]
    }
   ],
   "source": [
    "Hs=10 #cutting height of stacks\n",
    "#Input the TEDS csv file\n",
    "try:\n",
    "  df = read_csv('point.csv', encoding='big5')\n",
    "except:\n",
    "  df = read_csv('point.csv')\n",
    "# check_NOPandSCC(0)\n",
    "df = check_nan(df)\n",
    "# check and correct the X coordinates for isolated islands\n",
    "df = check_landsea(df)\n",
    "df = Elev_YPM(df)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58b298cd",
   "metadata": {},
   "source": [
    "- 使用`Hs`進行篩選「高空」點源"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "ac20b3e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "df=df.loc[(df.HEI>=Hs) & (df.NO_S.map(lambda x:x[0]=='P'))].reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcec9aef",
   "metadata": {},
   "source": [
    "- 排除沒有SNCPV排放者"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b93a2714",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['SUM']=[i+j+k+l+m for i,j,k,l,m in zip(df.SOX_EMI,df.NOX_EMI,df.CO_EMI,df.PM_EMI,df.NMHC_EMI)]\n",
    "df=df.loc[df.SUM>0].reset_index(drop=True)\n",
    "df['CP_NO'] = [i + j for i, j in zip(list(df['C_NO']), list(df['NO_S']))]\n",
    "df['DY1']=[i*j for i,j in zip(df.DW1,df.WY1)]\n",
    "df['HY1']=[i*j for i,j in zip(df.HD1,df.DY1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0f5922d",
   "metadata": {},
   "source": [
    "### 讀取並填滿CEMS資料檔\n",
    "- 填滿資料表\n",
    "  - 程式運作需要每筆**管煙**(**管編**+**煙編**)、每個小時都要有數值。將DataFrame轉成矩陣，再轉回DataFrame即可。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a5563fd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#71 factories with CEMS will emit (at ground) when stacks are operating\n",
    "fname=P+'point_cems.csv'\n",
    "cems=read_csv(fname)\n",
    "val='SOX PM NOX FLOW X_BLANK1 X_BLANK2'.split()\n",
    "nval=len(val)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "750bbd1d",
   "metadata": {},
   "source": [
    "- **管煙**欄位為處理過後檔案的特徵\n",
    "  - 新增**管編**+**煙編**(**管煙**)之新標籤、新增PM(設為SN之平均值)\n",
    "  - 新增時間標籤`MDH`。有的年度提供的CEMS檔案小時標記為0000\\~2300，有的是00\\~23，因此需要判斷一下\n",
    "  - 一個(**管煙**,**時籤**)組合只應對應一組CEMS數據，以`pivot_table sum`進行整併  \n",
    "- 維度標籤之計算。如果使用`標籤=序列.index(值)`指令將會非常耗時，直接使用`dict{值:標籤}`，會快很多。\n",
    "- `DataFrame`轉成`Array`，因所有值預設為0，如此就補滿空缺值了。有數據的部分再填入`Array`。\n",
    "- Array再轉回`DataFrame`。因後面CEMS數據使用邏輯太多樣了，`Array`形式不敷應用。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "dc63bb42",
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'CP_NO' not in cems.columns: #pre-process\n",
    "  cems=cems.drop(cems.loc[cems.C_NO=='C_NO'].index).reset_index(drop=True)\n",
    "#- 新增**管編**+**煙編**(**管煙**)之新標籤、新增PM(設為SN之平均值)\n",
    "\n",
    "  cems['CP_NO'] = [i + j for i, j in zip(list(cems['C_NO']), list(cems['NO_S']))]\n",
    "  cems['PM']=[(i+j)/2 for i,j in zip(cems.SOX,cems.NOX)]\n",
    "#- 新增時間標籤`MDH`。有的年度提供的CEMS檔案小時標記為0000~2300，有的是0~23，因此需要判斷一下\n",
    "\n",
    "  if max(cems.HOUR)>100:\n",
    "    cems['MDH']=[int(i*10000+j*100+k/100) for i,j,k in zip(cems.MONTH,cems.DATE,cems.HOUR)]\n",
    "  else:\n",
    "    cems['MDH']=[int(i*10000+j*100+k) for i,j,k in zip(cems.MONTH,cems.DATE,cems.HOUR)]\n",
    "#- 一個(**管煙**,**時籤**)組合只應對應一組CEMS數據，以`pivot_table sum`進行整併  \n",
    "\n",
    "  cems=pivot_table(cems,index=['CP_NO','MDH'],values=val,aggfunc=sum).reset_index()\n",
    "#- 維度標籤之計算。如果使用`標籤=序列.index(值)`指令將會非常耗時，直接使用`dict{值:標籤}`，會快很多。\n",
    "\n",
    "  #cems(df) convert to cemsM(matrix)\n",
    "  for MC in ['CP_NO','MDH']:\n",
    "    mc=MC.lower()\n",
    "    exec(mc+'=list(set(cems.'+MC+'))');exec(mc+'.sort()')\n",
    "    exec('n'+MC+'=len('+mc+')')\n",
    "    exec('d'+MC+'={'+mc+'[i]:i for i in range(n'+MC+')}')\n",
    "    exec('cems[\"i'+MC+'\"]=[d'+MC+'[i] for i in cems.'+MC+']')\n",
    "  if len(mdh)!=ndays*24:sys.exit('mdh coverage not enough!')    \n",
    "#- `DataFrame`轉成`Array`，因所有值預設為0，如此就補滿空缺值了。有數據的部分再填入`Array`。\n",
    "\n",
    "  cemsM=np.zeros(shape=(nMDH,nCP_NO,nval))\n",
    "  for i in range(nval):\n",
    "    cemsM[cems.iMDH[:],cems.iCP_NO[:],i]=cems[val[i]]\n",
    "#- Array再轉回`DataFrame`。因後面CEMS數據使用邏輯太多樣了，`Array`形式不敷應用。\n",
    "\n",
    "  DD={}\n",
    "  for i in range(nval):\n",
    "    DD[val[i]]=cemsM[:,:,i].flatten()\n",
    "  DD['MDH']  =[i for i in mdh for j in cp_no]\n",
    "  DD['CP_NO']=[j for i in mdh for j in cp_no]\n",
    "  cems=DataFrame(DD)\n",
    "  cems['C_NO']=[i[:8] for i in cems.CP_NO]\n",
    "  cems['MD']=[i//100 for i in cems.MDH]\n",
    "  cems.set_index('CP_NO').to_csv(fname)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f7a116b",
   "metadata": {},
   "source": [
    "### 整理CEMS各廠運作時間模式\n",
    "- 讀取資料庫欄位(**管編**`C_NO`, **管煙**`CP_NO`, 月日時`MDH`, 月日`MD`,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "80de4feb",
   "metadata": {},
   "outputs": [],
   "source": [
    "for MC in ['CP_NO','MDH','MD','C_NO']:\n",
    "  mc=MC.lower()\n",
    "  exec(mc+'=list(set(cems.'+MC+'))');exec(mc+'.sort()')\n",
    "  exec('n'+MC+'=len('+mc+')')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10dacc4f",
   "metadata": {},
   "source": [
    "- 沒有CEMS之同家工廠其它煙道、或鄰近工廠，可能具有上下游連動關係，因此適用該CEMS的操作特性，此處先行備妥。\n",
    "  - 個別廠的日操作特性，分析全年SOX排放量之小時變化，按排放量排序，代表最有可能運作之小時。\n",
    "  - 先執行`pivot_table`加總"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "3adb8316",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Hour of Day pattern\n",
    "cems['HR']=[i%100 for i in cems.MDH]\n",
    "pv_cems1=pivot_table(cems,index=['C_NO','HR'],values='SOX',aggfunc=sum).reset_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0599189f",
   "metadata": {},
   "source": [
    "  - 以**管編**為索引之新的`DataFrame`，其內容為小時之序位"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "ef5db925",
   "metadata": {},
   "outputs": [],
   "source": [
    "cems_HROD=DataFrame({'C_NO':c_no})\n",
    "cems_HROD['SOX_HR_ODER']=0\n",
    "for ic in cems_HROD.index:\n",
    "  pv1=pv_cems1.loc[pv_cems1.C_NO==c_no[ic]]\n",
    "  pv3=pv1.sort_values('SOX',ascending=False).reset_index(drop=True)\n",
    "  cems_HROD.loc[ic,'SOX_HR_ODER']=''.join(['{:d} '.format(i) for i in pv3.HR])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c56388bc",
   "metadata": {},
   "source": [
    "- 同樣原理應用在全年的逐日變化，內容為序列`mdh`的標籤。排序依據全天的流量總計，值越高者表當天最可能運作。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "5e4acc8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#orders for DY1\n",
    "pv_cems2=pivot_table(cems,index=['C_NO','MD'],values='FLOW',aggfunc=sum).reset_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37a2554c",
   "metadata": {},
   "source": [
    "  - `MD`欄位由**月日**轉變成`mdh`的標籤，轉變需要2步驟，否則會出錯"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "cde8d01f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Indexing is an exhaustive process.\n",
    "iMD=[mdh.index(i*100) for i in pv_cems2.MD] #change the MMDD into index sequence among MMDD00's\n",
    "pv_cems2.MD=iMD"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a64aa54",
   "metadata": {},
   "source": [
    "  - 以**管編**為索引之新的`DataFrame`，其內容為**日期標籤之序位**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "34bf8ecb",
   "metadata": {},
   "outputs": [],
   "source": [
    "cems_DAOD=DataFrame({'C_NO':c_no})\n",
    "cems_DAOD['FLOW_DA_ODER']=0\n",
    "for ic in cems_DAOD.index:\n",
    "  pv1=pv_cems2.loc[pv_cems2.C_NO==c_no[ic]]\n",
    "  pv3=pv1.sort_values('FLOW',ascending=False).reset_index(drop=True)\n",
    "  cems_DAOD.loc[ic,'FLOW_DA_ODER']=''.join(['{:d} '.format(i) for i in pv3.MD])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "837d78da",
   "metadata": {},
   "source": [
    "### 套用CEMS或運作時間模式\n",
    "- 建立各**管編**座標值之資料表"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "04ee0786",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfxy=pivot_table(df,index='C_NO',values=['UTM_E','UTM_N'],aggfunc=np.mean).reset_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c8e4c8f",
   "metadata": {},
   "source": [
    "- `BLS`為污染項目與布林值的對照表。以直接提取資料庫中符合布林值、要處理的筆數。\n",
    "  - `{'NMHC':'PM', 'NMHC':'PM'}`因為`NMHC`及`CO`沒有`CEMS`數據，假設與`PM`、`NOX`一樣。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "db67c67f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#booleans for pollutant selection\n",
    "c2v={'NMHC':'PM','SOX':'SOX','NOX':'NOX','PM':'PM','CO':'NOX'} #point.csv vs cems.csv\n",
    "BLS={c:df[c+'_EMI']>0 for c in c2v}\n",
    "colT=['HD1','DY1','HY1']\n",
    "col=['C_NO','CP_NO','HD1','DY1','HY1']+[i for i in df.columns if 'EMI' in i]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f30056c",
   "metadata": {},
   "source": [
    "- `s`為物質種類，須由引數讀取，且限定在`BLS`的索引範圍(`c2v`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "bc4ce5a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#for spe in [s for s in [sys.argv[1]] if s in BLS]:\n",
    "for spe in ['NMHC']: # for eg.\n",
    "  dfV=df[col].loc[BLS[spe]].reset_index(drop=True)   "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74871d69",
   "metadata": {},
   "source": [
    "- 對全廠加總並形成新的資料庫`dfV`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "072de6f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "for spe in ['NMHC']: # for eg.\n",
    "  dfV1=pivot_table(dfV,index='CP_NO',values=spe+'_EMI',aggfunc=sum).reset_index()\n",
    "  dfV2=pivot_table(dfV,index='CP_NO',values=colT,aggfunc=np.mean).reset_index()\n",
    "  dfV=merge(dfV1,dfV2,on='CP_NO')\n",
    "  dfV['C_NO']=[i[:8] for i in dfV.CP_NO]\n",
    "  for c in colT:\n",
    "    dfV[c]=np.array(dfV[c],dtype=int)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1e1f8f6",
   "metadata": {},
   "source": [
    "- 對照資料庫(`a`)與cems(`b`)的**管編**`cp`\n",
    "  - `ab`為二者交集\n",
    "  - `b1`為`b-a`，有CEMS數據卻不在資料庫中、沒有座標\n",
    "  - `c1`為`ab-b1`(就等於`ab`)\n",
    "  - 列出所有`c1`的座標序列"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "63cd6baf",
   "metadata": {},
   "outputs": [],
   "source": [
    "for spe in ['NMHC']: # for eg.\n",
    "  a,b=list(set(dfV.C_NO)),list(set(cems.C_NO));a.sort();b.sort()\n",
    "  ab=[i for i in a if i in b]\n",
    "  cp=list(set(dfV.CP_NO))\n",
    "  cp.sort()\n",
    "  ons=np.zeros(shape=(len(cp),nMDH))#,dtype=int)\n",
    "  #other fatories without CEMS, take the nearest one\n",
    "  b1=set(b)-set(dfxy.C_NO) #cems factory but without UTM location\n",
    "  c1=[c for c in b if c not in b1 and c in a] #cems plant with X,Y\n",
    "  cemsX=np.array([list(dfxy.loc[dfxy.C_NO==c,'UTM_E'])[0] for c in c1])\n",
    "  cemsY=np.array([list(dfxy.loc[dfxy.C_NO==c,'UTM_N'])[0] for c in c1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50f06a71",
   "metadata": {},
   "source": [
    "- 逐一**管編**進行迴圈\n",
    "  - 如果工廠沒有CEMS(`c not in ab`)，則以最近的一家廠的日變化特徵為其**時變係數**，令`c_cems`為該廠**管編**\n",
    "  - 如果該廠**日期標籤之序位**不足全年日數(ndays)，也將補足缺漏日之標籤`list(s365-set(pv2MD))`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "efcfcf84",
   "metadata": {},
   "outputs": [],
   "source": [
    "for spe in ['NMHC']: # for eg.\n",
    "  #loop for every factories\n",
    "  \n",
    "    c_cems=c\n",
    "    if c not in ab:\n",
    "      x0,y0=list(dfxy.loc[dfxy.C_NO==c,'UTM_E'])[0],list(dfxy.loc[dfxy.C_NO==c,'UTM_N'])[0]\n",
    "      dist=(cemsX-x0)**2+(cemsY-y0)**2\n",
    "      idx=list(dist).index(min(dist))\n",
    "      c_cems=c1[idx]\n",
    "    pv2MD=np.array(list(cems_DAOD.loc[cems_DAOD.C_NO==c_cems,'FLOW_DA_ODER'])[0].split(),dtype=int)\n",
    "    if len(pv2MD)<ndays: pv2MD=np.array(list(pv2MD)+list(s365-set(pv2MD)))\n",
    "    df_cp=dfV.loc[dfV.C_NO==c].reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6458e94c",
   "metadata": {},
   "source": [
    "- 該廠所有**煙編**逐一進行  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "2019f851",
   "metadata": {},
   "outputs": [],
   "source": [
    "for spe in ['NMHC']: # for eg.\n",
    "  for c in [i for i in a if i not in b1][:1]: #for eg.\n",
    "    #loop for every NO_S in this factory\n",
    "    for p in set(df_cp.CP_NO):\n",
    "      ip=cp.index(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdd11e10",
   "metadata": {},
   "source": [
    "  - 如果是CEMS煙道，其全年變化即為CEMS資料表中數值\n",
    "  - **point_cems.csv**檔案內容須事先確認其單位，此處要求其全年總合為1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "91f55c01",
   "metadata": {},
   "outputs": [],
   "source": [
    "for spe in ['NMHC']: # for eg.\n",
    "  for c in [i for i in a if i not in b1][:1]: #for eg.\n",
    "    #loop for every NO_S in this factory\n",
    "    for p in set(df_cp.CP_NO):\n",
    "      if p in set(cems.CP_NO):\n",
    "        ons[ip,:]=cems.loc[cems.CP_NO==p,c2v[spe]]*nhrs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e00b9db1",
   "metadata": {},
   "source": [
    "  - 將工廠操作矩陣(維度為`(日數、時數)`)先設定為0\n",
    "  - 如果不是，則由資料庫中讀取工作日數及小時數\n",
    "    - 如果是全天連續操作，則逐時都需標籤"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "7e99c569",
   "metadata": {},
   "outputs": [],
   "source": [
    "for spe in ['NMHC']: # for eg.\n",
    "  for c in [i for i in a if i not in b1][:1]: #for eg.\n",
    "    #loop for every NO_S in this factory\n",
    "    for p in set(df_cp.CP_NO):\n",
    "      if not p in set(cems.CP_NO):# for eg\n",
    "        dy1=dfV.DY1[ip]\n",
    "        hd1=dfV.HD1[ip]\n",
    "        md3=pv2MD[:dy1]\n",
    "        days=np.zeros(shape=(dy1,hd1),dtype=int)\n",
    "        if hd1==24:\n",
    "          hrs=np.array([i for i in range(24)],dtype=int)        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d12bdf9",
   "metadata": {},
   "source": [
    "  - 不是24小時，最有可能的小時，則由前面準備好的`cems_HROD.SOX_HR_ODER`第1小時開始，依序填入工作時數。\n",
    "    - 工作日則由為前述最可能工作的**日期標籤之序位**依序填入，形成`days`(月、日、時標籤)\n",
    "    - 將`days`壓平成為一維矩陣\n",
    "    - 這些有運作的時間，其`ons`值為1(其餘內設為0)   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "663e49dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "for spe in ['NMHC']: # for eg.\n",
    "  for c in [i for i in a if i not in b1][:1]: #for eg.\n",
    "    #loop for every NO_S in this factory\n",
    "    for p in set(df_cp.CP_NO):\n",
    "      if not p in set(cems.CP_NO):# for eg\n",
    "        if hd1!=24: #for eg.\n",
    "          first=np.array(list(cems_HROD.loc[cems_HROD.C_NO==c_cems,'SOX_HR_ODER'])[0].split(),dtype=int)[0]\n",
    "          hrs=np.array([(first+ih)%24 for ih in range(hd1)])\n",
    "        for id in range(dy1):\n",
    "          days[id,:]=md3[id]+hrs[:]\n",
    "        idx=days.flatten()\n",
    "        ons[ip,idx]=1.            "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3ff2cf0",
   "metadata": {},
   "source": [
    "### 輸出結果 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c3aa3f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "  fnameO=spe+'_ECP'+str(len(cp))+'_MDH'+str(len(mdh))+'_ONS.bin'\n",
    "  with FortranFile(fnameO, 'w') as f:\n",
    "    f.write_record(cp)\n",
    "    f.write_record(mdh)\n",
    "    f.write_record(ons)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}