鳥哥私房菜 - AERMOD - AERMAP 地形與受體

最近更新時間： 2019/08/07

如同前面的章節介紹的，整個擴散模式最重要的部份就是氣象資料。除了氣象資料之外，由於擴散還是有空間概念的，地表並不是平整的。比較高且接近煙流中心點的高度，會有比較高的濃度。因此，地表的著地濃度也需要考量高度的效應。此外，地球是圓的，因此表面並不是平面喔！但是，空氣品質模擬的時候，我們會假設地表示平整的。那如何處理地表的座標呢？這時就得要考量是需要使用經緯度座標？還是需要使用 UTM 座標？這都是好有趣的課題喔！

地理座標系與座標轉換
台灣地區的地表高程
AERMAP 的編譯與執行
參考資料

地理座標系與座標轉換

我們在數學上課程上面經常會有所謂的『座標』資料，包括平面座標、三維直角座標、極座標與球座標等等 (ps1)。不過，在實際應用中，大概使用到的通常是所謂的球體座標，跟經緯度有關的座標資料！這是因為我們目前經常使用的服務中，例如 google 地圖，就是使用經緯度座標來找出定位，讓我們可以經過 GPS 系統來到達目的地這一目標的緣故。不過，其實我們在短距離的環境中，比較習慣使用的還是屬於三維的直角座標系統啦！就是一般水平的 X, Y 軸，搭配高度的 Z 軸，配合成為直角座標的意思。直角座標系就有點像底下的圖示，包括平面與空間的直角座標。

地理座標系統

不過，地球是橢圓型球體，以整個地球的角度來看，似乎很難作為直角座標系統來處理！為了解決這個問題，因此，後來發展出經緯度座標(ps2)。以地球北極到南極畫出表面曲線的經度來說，0 度線是經過英國倫敦格林威治天文台舊址的線段，面向北極的右手邊，就稱為東經方向，左手邊則是西經方向。因為圓形是 360 度角，因此，向東或向西走 180 度就會碰面啊！因此，東經 180 度就是西經 180 度。至於緯度，則由赤道作為 0 度，向北為北緯，向南為南緯，緯度是地球某一點與地心連線對赤道的夾角所算出來的。在紀錄經緯度數值時，東經為 E，西經為 W，北緯為 N，南緯為 S。台灣地區大概位於東經 121 度以及北緯 25 度左右。

不過，這種經緯度座標的展示結果，可能會有很嚴重的誤差！舉例來說，如下圖左側，那個是將經緯度硬板成直角座標，從全球經緯度座標的直角座標中，你會發現，整個西伯利亞以及北極海都相當相當的大！南極冰層也很大！但如果你用球體展示，如下圖右方，就可以發現，其實北極冰層並沒有左側顯示的那樣龐大！這就是經緯度座標的問題！也就是說，同樣都是東經 120 到 121 度之間，但是在緯度 25 度以及緯度 50 度，這之間的距離就是不一樣長！因為不是直角座標啦！所以，看實際的距離，不能直接從經緯度去抓啦！可能要經過轉換成為距離單位 (公里) 之類的，會比較適當！

根據經緯度所設計出來的一套地理資訊系統，被稱為世界大地測量系統 (World Geoditic System, WGS)，這套系統主要用於地圖學、大地測量學與導航系統當中，從過去的 WGS60, WGS66, WGS72，到現在最新的 WGS84，都是透過經緯度座標來進行定位的一套系統。目前大部分的導航與定位系統，都是參考 WGS84 的參考座標來處理的。

UTM 座標 (通用橫軸墨卡托投影, Universal Transverse Mercator coordinate system)

如同前面提到的大地測量座標，大部份是以經緯度來處理。不過，經緯度在不同的國家地區，就是會有不一樣的長度情形！因此，就有許多的投影法再將地表處理成為類似直角座標的系統。其中一個很通用的座標系統，就是 UTM 座標了。

UTM 座標(ps4)以類似圓柱體包著地球，並將在圓柱體上小範圍的垂直線壓在地球表面上，然後將這個區塊再細分成更小的範圍，之後假設這範圍內的方塊形狀不變，來進行地圖的設計與繪製。為了將地球分為各個區間，所以，從西經 180 度的位置開始向東設計，每 6 度經度設計為一個 UTM 區域，所以共設計了 60 個 UTM 的區域。然後根據不同的緯度再設計小範圍的位置，從南緯 80 度到北緯 84 度之間，每 8 度給予一個區域，並加以編號 (從 C 開始)，就變成了一格一格的方塊格子了。

以台灣所在的位置來說，大概就是位於 50, 51 區的 Q, R 方塊左右。不過，傷腦筋的是，台灣地區剛剛好就在四個分區 (50Q, 51Q, 50R, 51R) 的正中央附近，這樣就無法使用正常的 UTM 座標分區來處理座標的運算了。那麼辦？後來在 1974 年左右，政府要繪製台灣的基本圖的時候，就以東經 120 ~ 122 橫跨兩個經度的範圍去設計 UTM 分區，這就被稱為是『二度分帶』的由來！只是，如此一來，台灣本島與澎湖就分屬不同的區了！所以，我們在參考圖示時，經常發現澎湖、台灣的 UTM 地圖是分開的！大概就是這麼回事！

但是，如果使用到這個區域，並且完全參考 UTM 原始座標的話，那麼我們台灣地區的許多位置，很可能就會變成負值。為了讓所有的座標都是正值，因此上述台灣的二度分帶 UTM 座標，就從原本的 51 區的原點向左移動經度一度，大約是 250 公里 (250000 公尺) 的距離，所以，原點的 X 軸就大概位於經度 120 度上面，那原點的 Y 軸使用赤道的位置，所以 UTM 座標的數值，就這樣可以算出來了。此外，UTM 座標的數值，大部分使用的是公尺為單位，你也可以使用公里為單位去撰寫 UTM 座標！根據這個座標，我們就可以使用類似直角座標系統來處理相關位置的計算了！

EPSG 與 Proj. 4

如上所述，地理資訊系統實在很複雜，有 WGS72, WGS84, UTM 座標，台灣地區甚至還有 TWD67 與 TWD97 的差異等等，如果沒有一個轉換機制，恐怕未來在互相參照彼此的座標時，會發生比較嚴重的問題。為此，EPSG 資料庫 (ps6) 就應運而生。 EPSG 官網所提供的地理資訊轉換系統大致上有微軟的 Access 、MySQL、PSQL 等資料庫的支援，有興趣的話，大家可以自行去下載實驗看看。 EPSG 針對不同的座標系統進行一些基礎的編號定義，我們以澎湖以及台灣的 TWD97 為例，這兩個 UTM 分區分別位於 (ps7)：

EPSG:3825 --> 澎湖 TM2, TWD97 二度分帶 (X 原點在東經 118 度)
EPSG:3826 --> 台灣本島 TM2, TWD97 二度分帶 (X 原點在東經 120 度)
EPSG:4326 --> WGS84 經緯度
EPSG:3824 --> 原本的 TWD97 經緯度

那如何進行座標轉換呢？除了 EPSG 已經具有某些轉換座標的特性之外，我們還能夠透過 PROJ (ps8) 這個軟體來進行座標軸的轉換喔！只是，這個 proj 軟體得要透過 EPEL 額外的軟體庫來安裝～因此，你需要在你的 CentOS Linux 上面安裝好軟體之後，才能夠進行座標軸的轉換。

# 1. 先使用 root 的身份，進行 epel 的安裝與 proj 的安裝
[root@localhost ~]# yum install epel-*
[root@localhost ~]# yum --enablerepo=epel install proj proj-epsg

# 2. 使用任何身份都可以，透過 cs2cs (coordinate system to coordinate system) 進行轉換：
[user@localhost ~]$ echo "121 25" | cs2cs +init=epsg:4326 +to +init=epsg:3826
250000.00       2765777.56 0.00

[user@localhost ~]$ echo "120.216 22.999" | cs2cs +init=epsg:4326 +to +init=epsg:3826
169628.02       2544387.25 0.00
# 這個座標點大概位於成大附近，這樣就可以進行經緯度轉 UTM 座標了！

先安裝好 proj 以及 proj-epsg 這兩個軟體，這樣就可以處理透過標準 EPSG 的資料庫內容轉換功能，可以直接轉換座標！處理的方式，就是先使用 echo 這個指令將你的經緯度座標 X 與 Y 整理好，然後帶入 cs2cs 這個指令功能即可！那如果想要從 UTM 座標轉成經緯度呢？就反過來寫即可！

# 1. 使用預設的格式，將 UTM 座標轉換成為經緯度座標：
[user@localhost ~]$ echo "169628 2544387" | cs2cs +init=epsg:3826 +to +init=epsg:4326
120d12'57.599"E 22d59'56.392"N 0.000
# 因為經緯度座標主要是以 60 度為分析，因此顯示的是 60 度角度的方式！

# 2. 使用如同前一個範例使用的經緯度座標，不要使用 60 度角度的方式處理：
[user@localhost ~]$ echo "169628 2544387" | cs2cs -f '%.4f' +init=epsg:3826 +to +init=epsg:4326
120.2160        22.9990 0.0000
# 增加了 -f '%.4f' ，代表使用小數點下 4 位數顯示的方式來處理的意思

透過這個 proj 的軟體，並透過預設的 EPSG 座標轉換標準，就可以順利的將台灣地區的 UTM 座標與經緯度座標進行切換了！這樣，未來我們在進行空品模式的模擬時，也比較容易判斷地理位置圖喔！

台灣地區的地表高程

因為地表高程會影響到模式模擬，過去台灣經歷多次的地表高程製圖，不過大部份都是空拍圖，實用性比較沒有這麼好。近年來因為 GPS 以及相關 GIS 等軟硬體較為發達，因此台灣內政部於 2016 年提供了全台灣 20 公尺一個格點的高度資料，使用了數值地形模型 (DTM) 的數值型態，提供給大眾使用，你可以從政府資料開放平台，或內政部地理資訊服務平台直接下載 (ps9)。

鳥哥是直接上內政部網站去下載的，如下連結所示。連接上之後，點選『2016年全臺灣20公尺網格DTM資料』的超連結，會出現另外一個視窗，鳥哥需要全台灣的資訊，因此鳥哥點選的是『不分幅全台及澎湖』那一橫列的『下載』按鈕。點選之後，瀏覽器就會跳出一個讓你選擇檔名的視窗。鳥哥輸入的是 taiwan.zip 這樣的檔名。最終該檔案就會被下載成為你剛剛輸入的檔名。請將這個檔案 (zip 壓縮格式) 傳送到你的伺服器上面去。

https://www.tgos.tw/TGOS/Web/MetaData/TGOS_Query_MetaData.aspx?key=TW-06-301000000A-612640

地形高度檔的解壓縮與轉檔格式

現在，假設我們將剛剛的 taiwan.zip 檔案放置到家目錄，那因為我們想要做的其實是 AERMOD 的附屬模式，也就是 AERMAP 這個模組。但是內政部提供的 TIFF 格式檔，但是 AERMAP 能讀取的則是需要 DEM 格式。也就是說，內政部提供的原始數據資料，並沒有辦法直接交給 AERMAP 去讀取的。那該如何是好？(ps10)

自由軟體 GDAL 提供了一個解決方案，那就是讓各種不同的地形資料進行資料的轉換！不過，根據測試，直接由 tiff 格式轉成 dem 格式，似乎會出問題。這是因為內政部提供的 tiff 檔案格式中，使用的是 EPSG:3826 的標準 TM2 TWD97 的格式，不過，許多方案使用的座標轉換，似乎都是僅支援標準的 UTM 60 zone 的分區 (還記得上面剛剛講過的 UTM 共 60 個區域的概念？)。另外，等等我們要使用的 AERMAP 軟體，似乎也是僅支援 60 個 UTM 分區而已。所以，使用標準正確的 EPSG:3826 反而無法轉換成功。

為此，我們就『假裝』內政部提供的座標軸資料是標準的 UTM 51 zone 這個區域 (其實台灣橫跨 50, 51 兩區喔)，因此預計使用 EPSG:32651 這個標準分區。基本上，北緯地區的座標參照，就是以 EPSG:326XX 來作為 UTM 座標的轉換，因為我們預計使用的是 zone 51，所以就將 XX 帶入 51 ，因此就得到 EPSG:32651 這個標準代碼了！詳細的各個代碼資訊，在你安裝過 proj-epsg 軟體後，前往觀察 /usr/share/proj/epsg 這個檔案，就可以得到正確的說明了。

我們預計要使用 GDAL 軟體來進行轉換，你得要使用 root 的身份來安裝好這個軟體。安裝完畢後，到家目錄去建立名為 taiwan_map 的目錄，並將剛剛從內政部的壓縮檔在該目錄下解開。最後進行轉檔，其中轉檔使用的表頭說明為 UTM 51 zone 的定義 (假的！為了符合給 AERMAP 使用！)。現在就一步一步慢慢來進行吧！

# 1. 先以 root 的身份來安裝好 gdal 軟體才行！
[root@localhost ~]# yum --enablerepo=epel install gdal
[root@localhost ~]# gdal_translate --help
Usage: gdal_translate [--help-general] [--long-usage]
       [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
             CInt16/CInt32/CFloat32/CFloat64}] [-strict]
       [-of format] [-b band] [-mask band] [-expand {gray|rgb|rgba}]
       [-outsize xsize[%] ysize[%]]
       [-unscale] [-scale[_bn] [src_min src_max [dst_min dst_max]]]* [-exponent[_bn] exp_val]*
       [-srcwin xoff yoff xsize ysize] [-projwin ulx uly lrx lry] [-epo] [-eco]
       [-a_srs srs_def] [-a_ullr ulx uly lrx lry] [-a_nodata value]
       [-gcp pixel line easting northing [elevation]]*
       [-mo "META-TAG=VALUE"]* [-q] [-sds]
       [-co "NAME=VALUE"]* [-stats] [-norat]
       src_dataset dst_dataset
# 出現這些資訊，代表 gdal 軟體已經安裝妥當了！

# 2. 使用一般帳號的身份，開始到 ~/taiwan_map 目錄底下去處理轉檔的事情。
[user@localhost ~]$ mkdir taiwan_map
[user@localhost ~]$ cd taiwan_map
[user@localhost taiwan_map]$ mv ~/taiwan.zip .  <==鳥哥假設，你的 taiwan.zip 放在你的家目錄底下喔！
[user@localhost taiwan_map]$ ll
-rw-rw-r--. 1 user user 104264769  8月  3 15:03 taiwan.zip  <==從內政部網站下載的全台地形檔

[user@localhost taiwan_map]$ unzip taiwan.zip
[user@localhost taiwan_map]$ ll
-rw-rw-r--. 1 user user        92  5月  9  2014 dem_20m.tfw
-rw-rw-r--. 1 user user 389235991 11月  2  2016 dem_20m.tif     <==台灣的地形高度
-rw-rw-r--. 1 user user        91  7月 24  2018 manifest.csv
-rw-rw-r--. 1 user user      8826  5月 31  2018 Metadata.xml
-rw-rw-r--. 1 user user        92  5月 22  2014 Penghu_20m.tfw
-rw-rw-r--. 1 user user   1045916  7月  7  2014 Penghu_20m.tif  <==澎湖的地形高度
-rw-rw-r--. 1 user user 104264769  8月  3 15:03 taiwan.zip

# 3. 開始以 EPSG:32651 進行轉檔的座標參考，記得要使用 USGSDEM 格式來轉檔喔！
[user@localhost taiwan_map]$ gdal_translate -a_srs EPSG:32651 -of USGSDEM dem_20m.tif taiwan_utm51.dem
Input file size is 10175, 19112
# 因為地形檔比較大，轉檔會花去幾分鐘的時間！最後會產生一個 1.1 G 左右的檔案！

[user@localhost taiwan_map]$ ll
-rw-rw-r--. 1 user user         92  5月  9  2014 dem_20m.tfw
-rw-rw-r--. 1 user user  389235991 11月  2  2016 dem_20m.tif
-rw-rw-r--. 1 user user         91  7月 24  2018 manifest.csv
-rw-rw-r--. 1 user user       8826  5月 31  2018 Metadata.xml
-rw-rw-r--. 1 user user         92  5月 22  2014 Penghu_20m.tfw
-rw-rw-r--. 1 user user    1045916  7月  7  2014 Penghu_20m.tif
-rw-rw-r--. 1 user user        388  8月  3 21:29 taiwan.prj
-rw-rw-r--. 1 user user 1177370624  8月  4 14:30 taiwan_utm51.dem            <==就是這個檔案囉！
-rw-rw-r--. 1 user user        177  8月  4 14:30 taiwan_utm51.dem.aux.xml
-rw-rw-r--. 1 user user  104264769  8月  3 15:03 taiwan.zip

如上所述，我們就會得到一個名為 taiwan_utm51.dem 的地形檔，要注意的是，其實內部的座標是使用 TWD97 的 UTM 格式，也就是 EPSG:3826 的標準座標。不過因為要給 AERMAP 使用，所以，表頭資料給的則是 EPSG:32651 的 UTM 51 zone。事實上，內部的座標依舊採用 EPSG:3826 的數值！由於與標準的定義不同，所以這裡得要持續特別強調喔！

台灣地圖使用範圍

因為 TM2 TWD97 的區域是二度分帶，台灣本島地區位於東經 120 ~ 122 度之間，超過這個範圍就無法處理了。同時，因為地圖僅有台灣地區 (上述這個 DEM 檔案內容)，經過分析後，這個地形高度檔案的 UTM 範圍大致上是這樣的 (單位為公尺)：

X 軸的範圍： 148320 ~ 351800
Y 軸的範圍： 2419500 ~ 2801720

當你設計模擬範圍的時候，不要超出這個範圍之外喔！否則 AERMAP 可能會無法運作的！

AERMAP 的編譯與執行

如前所述，AERMOD 要執行的時候，還得要搭配地表高程的 AERMAP 模組去處理出所需要的模擬範圍區間的地形。而上面我們已經做好了 AERMAP 能夠處理的『假的 UTM 51 zone 的台灣 UTM 座標』資料，接下來就來處理 AERMAP 的相關處理流程！第一個需要處理的，當然就是 AERMAP 的編譯。

編譯 AERMAP

AERMAP 這次鳥哥也是使用 v18081 版本，請自行下載，並在伺服器使用者家目錄的 aermap 目錄解開，之後再進行編譯處理的行為！整個流程比 AERMET 要簡單！因為大部分的編譯程式，已經在前面安裝過，所以，這裡就簡單處理即可：

# 1. 在家目錄底下建立 aermap 目錄，進入該目錄後，開始下載 aermap 程式碼
[user@localhost ~]$ mkdir aermap
[user@localhost ~]$ cd aermap
[user@localhost aermap]$ wget https://www3.epa.gov/ttn/scram/models/aermod/aermap/aermap_source.zip
[user@localhost aermap]$ ll
-rw-rw-r--. 1 user user 199162  4月 24  2018 aermap_source.zip

[user@localhost aermap]$ mkdir build
[user@localhost aermap]$ cd build
[user@localhost build]$ unzip ../aermap_source.zip
[user@localhost build]$ ll
-rw-rw-r--. 1 user user 311998  4月 19  2018 aermap.f
-rw-rw-r--. 1 user user  33231  4月 23  2018 aermap_mcbs.txt
-rw-rw-r--. 1 user user   1647  4月 19  2018 gfortran-aermap-32bit.bat
-rw-rw-r--. 1 user user   1627  4月 19  2018 gfortran-aermap-64bit.bat  <==編譯參考檔
-rw-rw-r--. 1 user user   2467  4月 19  2018 intel-aermap.bat
-rw-rw-r--. 1 user user  32478  4月 19  2018 mod_main1.f
-rw-rw-r--. 1 user user   2713  4月 19  2018 mod_tifftags.f
-rw-rw-r--. 1 user user      8  2月 22  2018 README.md
-rw-rw-r--. 1 user user  49490  4月 19  2018 sub_calchc.f
-rw-rw-r--. 1 user user  30062  4月 19  2018 sub_chkadj.f
-rw-rw-r--. 1 user user  25862  4月 19  2018 sub_chkext.f
-rw-rw-r--. 1 user user  21866  4月 19  2018 sub_cnrcnv.f
-rw-rw-r--. 1 user user  28213  4月 19  2018 sub_demchk.f
-rw-rw-r--. 1 user user  37991  4月 19  2018 sub_demrec.f
-rw-rw-r--. 1 user user  38215  4月 19  2018 sub_demsrc.f
-rw-rw-r--. 1 user user  18813  4月 19  2018 sub_domcnv.f
-rw-rw-r--. 1 user user  28722  4月 19  2018 sub_initer_dem.f
-rw-rw-r--. 1 user user  49116  4月 19  2018 sub_initer_ned.f
-rw-rw-r--. 1 user user  43539  4月 19  2018 sub_nadcon.f
-rw-rw-r--. 1 user user  94251  4月 19  2018 sub_nedchk.f
-rw-rw-r--. 1 user user  54702  4月 19  2018 sub_read_tifftags.f
-rw-rw-r--. 1 user user  10504  4月 19  2018 sub_reccnv.f
-rw-rw-r--. 1 user user  44332  4月 19  2018 sub_recelv.f
-rw-rw-r--. 1 user user  10292  4月 19  2018 sub_srccnv.f
-rw-rw-r--. 1 user user  41882  4月 19  2018 sub_srcelv.f
-rw-rw-r--. 1 user user  15855  4月 19  2018 sub_utmgeo.f

[user@localhost build]$ vim compile.sh
#!/bin/bash

COMPILE_FLAGS=" -fbounds-check -Wuninitialized -Ofast -static -march=native -ffast-math -funroll-loops"
LINK_FLAGS=" -static -Ofast -march=native -ffast-math -funroll-loops"
sources="
 mod_main1.f
 mod_tifftags.f
 aermap.f
 sub_calchc.f
 sub_chkadj.f
 sub_chkext.f
 sub_demchk.f
 sub_nedchk.f
 sub_cnrcnv.f
 sub_demrec.f
 sub_demsrc.f
 sub_domcnv.f
 sub_initer_dem.f
 sub_initer_ned.f
 sub_nadcon.f
 sub_reccnv.f
 sub_recelv.f
 sub_srccnv.f
 sub_srcelv.f
 sub_utmgeo.f
 sub_read_tifftags.f"

objects=$( echo ${sources} | sed 's/\.f/\.o/g' )

for file in ${sources}
do
        echo "gfortran -m64 -c ${COMPILE_FLAGS} ${file}"
        gfortran -m64 -c ${COMPILE_FLAGS} ${file}
done

echo "gfortran -m64 -o aermap.exe ${LINK_FLAGS} ${objects}"
gfortran -m64 -o aermap.exe ${LINK_FLAGS} ${objects}

[user@localhost build]$ sh compile.sh
[user@localhost build]$ ll *.exe
-rwxrwxr-x. 1 user user 2887656  8月  4 21:20 aermap.exe
# 這樣就編譯成功了！這個程式準備開始使用囉！

AERMAP 概說

AERMAP 主要是用來作為 AERMOD 的地形資料附屬模式，主要的目的就是取得 AERMOD 運算時所需要的地形資訊。那 AERMAP 必須要有實際測量而來的地形資料才行啊！AERMAP 預設讀取的兩種格式有：

DEM 格式：主要是以 USGS DEM 格式儲存的地表高度資訊 (就是我們之前轉換的格式)
美國 NED 資料，使用 GeoTIFF 格式。但只限於美國的 NED 資料，因此，我們台灣地區當然忽略囉！

除了特殊的 NED 格式的資料之外，AERMAP 主要支援的座標系統為 UTM 座標，如前所述，整個地球共區分為 60 個 UTM 區域，台灣地區位於 50 到 51 之間。但台灣地區所使用的座標系統並不在標準的 60 個 UTM 裡面，而是使用 TM2 TWD97 的 UTM 座標，是 2 度分帶座標，而不是傳統的 6 度分帶 UTM 座標。且原點並不是在原有的 UTM 標準座標內。但由於 AERMAP 好像僅支援標準 UTM 分區，所以，我們才會在前面的小節，假裝台灣的 TM2 TWD97 座標為 UTM 51 區！事實上還是使用標準 TM2 TWD97 的座標資料啦。

基本上，AERMAP 除了 UTM 座標之外，也是支援經緯度座標的，但這當然得要視你的輸入資料而定。那如何指定是 UTM 座標還是經緯度座標呢？在 AERMAP 的設定檔裡面，可以使用 DOMAINXY (UTM 的 X, Y 座標) 或 DOMAINLL (經緯度的 Lat, Lon 座標) 來規範。

要注意的是，越多的受體點，就需要越多的 AERMOD 運算時間！所以，雖然我們使用的台灣地區地形高程資料是細到 20 公尺的平均高度，但是，如果要模擬的區域全部 50km x 50km 的區域全部都以 20m 來作為一個格點計算，那就會產生至少六百萬個受體點！光是 AERMOD 模式輸出的資料，恐怕就要讓我們等到頭昏腦脹～因此，最好根據你的煙囪位置，以及盛行風，還有你有興趣知道的受體點來規劃地形高度！這樣在最終輸出所有各個受體點的資料時，才不會這麼恐怖...而為了簡化受體點網格的設計，而有底下的幾種常見的網格設計方式：

如上所示，我們可以指定用放射狀的方式設計格點 (左下方)，也可以使用等距離的方式設計格點 (右下方)，更能夠使用單一格點來進行受體點的規範！這都是可以進行的。一般來說，在排放源中心可能會設計比較緊密的網格，或者是預期所在的著地點會進行較密集的格點設計，而在空地或比較不在乎的區域，可以將格點設計的寬鬆一些。

在單位的處理方面，雖然 AERMAP 可支援很多的長度單位，不過，最終輸出的單位是公尺喔！所以，建議所有的輸入單位都統一使用公尺，這樣在設計與分析數據上面，才不會突然間出問題。這一點對我們來說，還好，畢竟台灣地區通常使用的長度單位就是公尺或公里。這裡要統一喔～使用公尺！公尺！

另外，在設計 AERMAP 的設定檔時，幾個比較重要的項目如下所示。接下來我們分別介紹每個部份囉！

CO：就是 control 的意思，主要設計網格所在的 UTM 座標，以及是否需要調整座標格點的定位等
RE：就是 receptor 的意思，主要是用來設計受體點的座標位置！這就重要了！
OU：就是 output 的意思，設計輸出的檔名資料。

CO 參數，控制網格位置與 UTM 座標還有座標轉換

現在，假設我們需要模擬的是興達火力電廠的煙囪排放，那麼興達電廠的煙囪座標是哪裡？其實，除了查詢實際的排放量資料之外，我們可以透過 google map 的平面地圖與衛星地圖互相切換搜尋，可以找到煙囪大致上座標是在『 120.198721, 22.852061 』這個經緯度上面，透過 cs2cs 的轉換，這個座標大致的 UTM 座標是這樣的：

[user@localhost ~]$ echo "120.198721 22.852061" | cs2cs +init=epsg:4326 +to +init=epsg:3826
167767.93       2528124.87 0.00

座標大致上就是在 167768, 2528125 的 UTM 座標上！通常設計時，選個比較好處理的整數 (公里) 會比較容易撰寫設定檔，加上南台灣空氣品質不佳的時候，大部份是冬天，冬天吹東北季風，因此，我們將模擬區域的中心點，可以朝向東北方移動一點點，所以，預計使用的網格中心點設計為 168000, 2529000 。這是因為 UTM 座標軸跟我們一般認知的直角座標一樣，向右、向上為正，也就是向東、向北為正的意思。所以向東北方移動到最近的公里數值，就得到 168,000、2529,000 這個 X, Y 數值了。

中心點為： 168000.0 2529000.0、假設位於 UTM 51 zone 的位置上。
由於模擬 50 公里，所以由中心點向四個延伸 25 公里。由於模擬網格大多為長方形，因此只要紀錄左下角、右上角的座標即可。不過台灣地區地圖指到 148xxx 的網格點，所以這裡我使用了 149000 這個左方的格點喔！最終就可以得到左下、右上的座標分別是： (149000, 2504000), (193000, 2554000)。
由於需要套圖，所以使用的座標軸不需要轉換。

基本上，最陽春的 CO 參數有點像底下這樣：

[user@localhost ~]$ vim ~/aermap/aermap_xingda.inp
CO STARTING
   TITLEONE  terrain, 20m taiwan map in fake UTM 51 zone
   TITLETWO  this file calculate for xingda power plan.
   DATATYPE  DEM
   DATAFILE  taiwan_utm51.dem
   DOMAINXY  149000.0 2504000.0 51 193000.0 2554000.0 51
   ANCHORXY  168000.0 2529000.0 168000.0 2529000.0 51  0
   RUNORNOT  RUN
CO FINISHED

TITLEONE, TITLETWO：只是說明這個檔案的內容，也就是整體 AERMAP 模式的輸入、輸出資料說明而已。建議不要超過 68 個字元之外～手冊是這樣說的，可能有偵測的長度限制。在這個範例中，TITLEONE 我指定的是輸入資料的來源， TITLETWO 則是說明主網格在哪個地區。
DATATYPE：就是我們剛剛轉換出來的那個地形高度檔案，格式就是 DEM 喔！
DATAFILE：剛剛的檔名～我們會將這個檔案做個連結檔的格式放置到本目錄中
DOMAINXY：定義出左下角的 UTM X, Y 座標，以及 51 UTM 分區，以及右上角的 UTM X, Y 座標，還以及 51 UTM 51 分區。
ANCHORXY：是否需要進行座標軸轉換～如果不需要座標轉換，則輸入兩次中央的座標軸資料，然後輸入 51 0 (UTM 分區)，這樣就可以了。什麼時候需要座標轉換呢？如果你想要讓排放中心作為原點 (0, 0)，同時不用顯示地圖 (或者是地圖也已經經過座標轉換)，那就可以透過類似如下的方式來將中心點變成是原點的模樣：
```
   ANCHORXY  0.0 0.0 168000.0 2529000.0 51  0
```
那原本的 168000, 2529000 座標，就會偏移到變成 (0, 0) 這樣，所有的座標展示都會以 0, 0 作為偏移點！
RUNORNOT：設定是要跑模式還是不跑。當然就是 RUN 啊！

RE 參數，設計受體點的座標位置

關於受體點的設計，主要有平行網格點的設計，透過 GRIDCART 來處理，也可以透過放射狀展示的 GRIDPOLR 來處理，最後，也可以直接指定某個受體點，利用 DISCCART 來處理！相關的設計可以參考手冊 3.4.2 ~ 3.4.5 小節的說明，主要分佈在手冊 page 3-26 到 page 3-38 的內容來處理 (ps12)。鳥哥個人的設計比較習慣使用平行格點以及單一格點，所以，底下會以 GRIDCART 搭配 DISCCART 來說明喔！

至於格點的設計上，由於風向的轉變情況在全年可能會有南北對調的改變，因此，整個污染物可能會回吹到中心點附近。因此，越靠近中心，我們設計的格點會較細。現在假設距離中心點左右各 5 公里的網格，每 200 公尺一個格點，而在之外的地方，則是每 500 公尺一個格點。為了方便以矩形來設計格點，所以，我們以底下的圖示來設計所需要的受體點：

基本上共分為五個區域，只有 A 區是 200 公尺一個點，其他都是 500 公尺一個點。我們以左下角的點作為該區的原點，因此，這五個區域的特殊參數就會變成這樣(假設原本的中心點 UTM 座標為 Xori, Yori 喔！)：

A：X = (Xori - 5000), Y = (Yori - 5000), X 個數 (10000/200) = 50, Y 個數也是 50，格點距離 200 公尺
B：X = (Xori - 5000), Y = (Yori + 5000), X 個數 (10000/500) = 20, Y 個數 (20000/500) = 40，格點距離 500 公尺
C：X = (Xori - 5000), Y = (Yori - 25000), X 個數 (10000/500) = 20, Y 個數 (20000/500) = 40，格點距離 500 公尺
D：X = (Xori - 25000), Y = (Yori - 25000), X 個數 (20000/500) = 40, Y 個數 (50000/500) = 100，格點距離 500 公尺
E：X = (Xori + 5000), Y = (Yori - 25000), X 個數 (20000/500) = 40, Y 個數 (50000/500) = 100，格點距離 500 公尺

至於網格受體點的設計使用 GRIDCART 這個設定來處理，這個設定的參數有點像底下這樣：

   GRIDCART ID_name STA
                   XYINC Xinit Xnum Xdelta Yinit Ynum Ydelta
   GRIDCART ID_name END

ID_name：可以隨便你取，不過，習慣上，可能根據上面的設計，定義為 BLOCKA, BLOCKB, BLOCKC, BLOCKD 與 BLOCKE 這樣
Xinit：就是該網格點左下角的 UTM 座標 X 軸
Xnum：就是 X 軸上面的格點數量
Xdelta：就是每個格點的距離
Yinit, Ynum, Ydelta：跟上面的 X 類似，只是變成 Y 軸的資料而已。

不過，在我們的案例中，比較特別的是，因為興達的煙囪太左側了，所以左側的格點只到 168K-149K = 19K 而已，並不是 25 公里。另外，因為 UTM 能接受的領域範圍大概在 153K 以東，所以建議，D 區域的設計需要調整一下 (X 軸從 25K 變成 15K 的變化)。此外，我們並以中心點 (168000, 2529000) 設計出單一一個格點的中心位置，就透過 DISCCART 這個參數。所以，根據上面的說明，設計出來的 RE 參數，會有點像底下這樣：

[user@localhost ~]$ vim ~/aermap/aermap_xingda.inp
RE STARTING
   GRIDCART BLOCKA STA
                   XYINC 158000. 200 200. 2519000. 200 200.
   GRIDCART BLOCKA END
   GRIDCART BLOCKB STA
                   XYINC 158000. 40  500. 2529000.  30 500.
   GRIDCART BLOCKB END
   GRIDCART BLOCKC STA
                   XYINC 158000. 40  500. 2504000.  30 500.
   GRIDCART BLOCKC END
   GRIDCART BLOCKD STA
                   XYINC 149000. 18  500. 2504000. 100 500.
   GRIDCART BLOCKD END
   GRIDCART BLOCKE STA
                   XYINC 178000. 30  500. 2504000. 100 500.
   GRIDCART BLOCKE END
   DISCCART 168000.0 2529000.0
RE FINISHED

這樣就處理好我們所需要的網格資料，大概會計算出 10101 個網格點，數量非常龐大！因此，使用 AERMAP 去計算時，會花費一小段時間，大概幾分鐘吧！計算出來的格點，就會寫入到底下的 OU 參數設定當中囉！

OU 參數，只要設計好輸出檔案

最終，我們總是需要一個輸出檔，這個設計就很簡單：

[user@localhost ~]$ vim ~/aermap/aermap_xingda.inp
OU STARTING
   RECEPTOR  AERMAP_XINGDA.REC
OU FINISHED

習慣上面檔名還是會寫大寫，此外，副檔名鳥哥這邊寫了受體 (receptor) 的縮寫，就是 .REC 的模樣！接下來就開始執行他：

# 1. 先進入到 aermap 的目錄，然後將台灣的 DEM 格式地形檔做個連結過來：
[user@localhost ~]$ cd ~/aermap
[user@localhost aermap]$ ln -s ../taiwan_map/taiwan_utm51.dem .
[user@localhost aermap]$ ll
-rw-rw-r--. 1 user user     958  8月  7 12:08 aermap_xingda.inp
drwxrwxr-x. 2 user user    4096  8月  4 21:20 build
lrwxrwxrwx. 1 user user      30  8月  4 21:32 taiwan_utm51.dem -> ../taiwan_map/taiwan_utm51.dem

# 2. 開始執行模式的運作：
[user@localhost aermap]$ ./build/aermap.exe aermap_xingda.inp
 usage: 0, 1, or 2 args

 Usage: AERMAP 18081 takes either no or one or two parameters.
        Either
              AERMAP
        Or
              AERMAP plumetest.inp
        Or
              AERMAP plumetest.inp plumetest.out

        The first parameter  is the .INP file name,
        The second parameter is the .OUT file name,

+Now Processing SETUP Information
+Processing Setup...

***********************
 OPENING FILE:

   DEM File #:         1
   DEM File Name: taiwan_utm51.dem

   Partial analysis of file structure using first    10240 characters.

     Number of bytes read:        10240
     File Type:            One Continuous Record - ASCII File

 CLOSING THE FILE.

 Exiting DEMCHK
 Exiting CHKADJ
 Exiting DOMCNV
 Exiting RECCNV
 Exiting CHKEXT
 Exiting DEMREC

+Initializing Terrain Data...
 This may take few a minutes...

 Exiting INITER_DEM

+Now Processing Receptor       1 of   10101
+Now Processing Receptor       2 of   10101
+Now Processing Receptor       3 of   10101
....
+Now Processing Receptor   10100 of   10101
+Now Processing Receptor   10101 of   10101

[user@localhost aermap]$ ll
-rw-rw-r--. 1 user user    4105  8月  7 12:08 aermap_xingda_DOMDETAIL.OUT
-rw-rw-r--. 1 user user     958  8月  7 12:08 aermap_xingda.inp
-rw-rw-r--. 1 user user    3170  8月  7 12:08 aermap_xingda_MAPDETAIL.OUT
-rw-rw-r--. 1 user user    2448  8月  7 12:08 aermap_xingda_MAPPARAMS.OUT
-rw-rw-r--. 1 user user   91824  8月  7 12:12 aermap_xingda.out
-rw-rw-r--. 1 user user  299142  8月  7 12:12 AERMAP_XINGDA.REC   <==這就是輸出檔！
drwxrwxr-x. 2 user user    4096  8月  4 21:20 build
lrwxrwxrwx. 1 user user      30  8月  4 21:32 taiwan_utm51.dem -> ../taiwan_map/taiwan_utm51.dem

這個 AERMAP_XINGDA.REC 就是我們最終需要使用到的地形高度資料，只是目前暫時無法使用軟體繪製其高度而已。

參考資料

(ps1)座標系的數學式：https://zh.wikipedia.org/wiki/坐標系
(ps2)地理座標系統：https://zh.wikipedia.org/wiki/地理坐标系
(ps3)世界大地測量系統 (World Geodetic System, WGS)：https://zh.wikipedia.org/wiki/世界大地测量系统
(ps4)UTM座標的解釋：https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system
(ps5)上河文化『大地座標系統漫談』：http://www.sunriver.com.tw/grid_tm2.htm
(ps6)EPSG 官網：http://www.epsg.org/
(ps7)QGIS 資源網：http://gis.rchss.sinica.edu.tw/qgis/?p=2823
(ps8)使用 Proj 軟體 cs2cs 轉換座標：https://proj.org/apps/cs2cs.html
(ps9)台灣內政部提供的 20m 間隔的地形模型：https://data.gov.tw/dataset/35430
直接連結到內政部網站：https://www.tgos.tw/TGOS/Web/MetaData/TGOS_Query_MetaData.aspx?key=TW-06-301000000A-612640
(ps10)將各種不同的格式轉檔成為 DEM：https://library.carleton.ca/help/dem-formats#USGS_DEM
GDAL 轉換軟體：https://geekforum.wordpress.com/category/gdal/
(ps11)US EPA 的 AERMAP 模式相關資料：
https://www.epa.gov/scram/air-quality-dispersion-modeling-related-model-support-programs#aermap
(ps12)AERMAP 的使用者手冊：https://www3.epa.gov/ttn/scram/models/aermod/aermap/aermap_userguide_v18081.pdf
鳥哥這裡的備份：https://linux.vbird.org/enve/images/aermap_userguide_v18081.pdf