鸟哥私房菜 - 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
鸟哥这里的备份：http://vbird.org.cn/enve/images/aermap_userguide_v18081.pdf