空气质量仿真软件：CAMx (Comprehensive Air Quality Model with Extensions)_（9）.地形和网格设置

最新推荐文章于 2025-05-03 08:19:19 发布

kkchenjj

最新推荐文章于 2025-05-03 08:19:19 发布

阅读量733

点赞数 16

分类专栏：环境仿真文章标签：仿真模拟环境仿真模拟仿真

本文链接：https://blog.csdn.net/weixin_42749425/article/details/147675148

版权

环境仿真专栏收录该内容

255 篇文章

订阅专栏

地形和网格设置

地形数据的准备

在进行空气质量仿真时，地形数据是至关重要的输入之一。地形数据的准确性和分辨率直接影响到模型的模拟结果。CAMx 支持多种地形数据格式，常见的包括 DEM（Digital Elevation Model）和 TERR（Terrain Data）格式。地形数据通常需要从外部数据源获取，如 USGS（United States Geological Survey）或 GIS（Geographic Information System）软件。

在这里插入图片描述

1. 获取地形数据

1.1 从 USGS 获取 DEM 数据

USGS 提供了大量的地形数据，可以通过其官方网站下载。以下是一个简单的步骤说明：

访问 USGS 网站。
选择你感兴趣的区域。
下载所需的 DEM 数据文件，通常是 .hgt 或 .dem 格式。

1.2 从 GIS 软件获取 TERR 数据

使用 GIS 软件如 ArcGIS，可以生成 TERR 格式的地形数据。以下是一个简单的步骤说明：

打开 ArcGIS。
添加你要使用的 DEM 数据。
使用工具将 DEM 数据转换为 TERR 格式。例如，使用 Raster to ASCII 工具将 DEM 转换为 ASCII 格式，然后进一步转换为 TERR 格式。

2. 地形数据的处理

2.1 使用 GDAL 工具处理地形数据

GDAL（Geospatial Data Abstraction Library）是一个强大的开源地理空间数据处理库，可以用于处理 DEM 数据。以下是一个示例代码，展示如何使用 GDAL 将 DEM 数据转换为 ASCII 格式：


# 导入 GDAL 库

from osgeo import gdal

import numpy as np



# 读取 DEM 数据

dem_file = 'path/to/your/dem/file.dem'

dataset = gdal.Open(dem_file)



# 获取地形数据

elevation_band = dataset.GetRasterBand(1)

elevation_data = elevation_band.ReadAsArray()



# 获取地理信息

geotransform = dataset.GetGeoTransform()

origin_x = geotransform[0]

origin_y = geotransform[3]

pixel_width = geotransform[1]

pixel_height = geotransform[5]



# 写入 ASCII 文件

output_file = 'path/to/output/terrain.asc'

header = f'ncols {elevation_data.shape[1]}\n'

header += f'nrows {elevation_data.shape[0]}\n'

header += f'xllcorner {origin_x}\n'

header += f'llcorner {origin_y}\n'

header += f'cellsize {pixel_width}\n'

header += 'NODATA_value -9999\n'



with open(output_file, 'w') as f:

    f.write(header)

    np.savetxt(f, elevation_data, fmt='%1.2f')

3. 网格设置

CAMx 中的网格设置是模型运行的基础。合理的网格设置可以提高模型的精度和计算效率。网格可以分为水平网格和垂直网格。

3.1 水平网格设置

水平网格设置通常包括网格分辨率、网格区域和网格投影。CAMx 支持多种投影方式，包括经纬度投影和 UTM 投影。

3.1.1 经纬度投影

经纬度投影是最常见的网格设置方式。以下是一个示例代码，展示如何在 CAMx 中设置经纬度投影的网格：


# 设置经纬度投影的网格

# 假设我们有一个区域，其范围为经度 110.0 到 120.0，纬度 30.0 到 40.0

# 水平网格分辨率为 0.1 度



# 定义网格参数

lon_min = 110.0

lon_max = 120.0

lat_min = 30.0

lat_max = 40.0

resolution = 0.1



# 计算网格尺寸

n_lon = int((lon_max - lon_min) / resolution) + 1

n_lat = int((lat_max - lat_min) / resolution) + 1



# 生成网格

lons = np.linspace(lon_min, lon_max, n_lon)

lats = np.linspace(lat_min, lat_max, n_lat)



# 保存网格参数到文件

grid_file = 'path/to/output/grid.txt'

with open(grid_file, 'w') as f:

    f.write(f'Longitude_min: {lon_min}\n')

    f.write(f'Longitude_max: {lon_max}\n')

    f.write(f'Latitude_min: {lat_min}\n')

    f.write(f'Latitude_max: {lat_max}\n')

    f.write(f'Resolution: {resolution}\n')

    f.write(f'Number_of_Longitudes: {n_lon}\n')

    f.write(f'Number_of_Latitudes: {n_lat}\n')

3.1.2 UTM 投影

UTM 投影（Universal Transverse Mercator）是一种常用的投影方式，适用于较小的区域。以下是一个示例代码，展示如何在 CAMx 中设置 UTM 投影的网格：


# 设置 UTM 投影的网格

# 假设我们有一个区域，其范围为 UTM 东经 500000 到 600000，北纬 3500000 到 3600000

# 水平网格分辨率为 1000 米



# 定义网格参数

utm_east_min = 500000

utm_east_max = 600000

utm_north_min = 3500000

utm_north_max = 3600000

resolution = 1000



# 计算网格尺寸

n_east = int((utm_east_max - utm_east_min) / resolution) + 1

n_north = int((utm_north_max - utm_north_min) / resolution) + 1



# 生成网格

eastings = np.linspace(utm_east_min, utm_east_max, n_east)

northings = np.linspace(utm_north_min, utm_north_max, n_north)



# 保存网格参数到文件

grid_file = 'path/to/output/utm_grid.txt'

with open(grid_file, 'w') as f:

    f.write(f'UTM_East_min: {utm_east_min}\n')

    f.write(f'UTM_East_max: {utm_east_max}\n')

    f.write(f'UTM_North_min: {utm_north_min}\n')

    f.write(f'UTM_North_max: {utm_north_max}\n')

    f.write(f'Resolution: {resolution}\n')

    f.write(f'Number_of_Eastings: {n_east}\n')

    f.write(f'Number_of_Northings: {n_north}\n')

3.2 垂直网格设置

垂直网格设置包括网格层数和层间高度。合理的垂直网格设置可以更好地模拟大气垂直结构。

3.2.1 垂直网格层数

CAMx 中的垂直网格层数可以根据研究目的进行设置。以下是一个示例代码，展示如何设置垂直网格层数：


# 设置垂直网格层数

# 假设我们设置 10 层垂直网格



# 定义垂直网格参数

n_layers = 10



# 生成垂直网格

layers = np.linspace(0, 1000, n_layers)



# 保存垂直网格参数到文件

vertical_grid_file = 'path/to/output/vertical_grid.txt'

with open(vertical_grid_file, 'w') as f:

    f.write(f'Number_of_Layers: {n_layers}\n')

    f.write('Layer_heights (m):\n')

    np.savetxt(f, layers, fmt='%1.2f')

3.2.2 垂直网格的高度分布

垂直网格的高度分布可以是均匀的，也可以是分层的。以下是一个示例代码，展示如何设置分层的垂直网格高度分布：


# 设置分层的垂直网格高度分布

# 假设我们设置 10 层垂直网格，高度分布为 [0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000] 米



# 定义垂直网格参数

n_layers = 10

layer_heights = [0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000]



# 保存垂直网格参数到文件

vertical_grid_file = 'path/to/output/vertical_grid.txt'

with open(vertical_grid_file, 'w') as f:

    f.write(f'Number_of_Layers: {n_layers}\n')

    f.write('Layer_heights (m):\n')

    for height in layer_heights:

        f.write(f'{height}\n')

4. 地形数据与网格的结合

将地形数据与网格相结合是空气质量仿真中的关键步骤。以下是一个示例代码，展示如何将 DEM 数据与水平网格相结合：


# 导入必要的库

import numpy as np

from osgeo import gdal



# 读取 DEM 数据

dem_file = 'path/to/your/dem/file.dem'

dataset = gdal.Open(dem_file)

elevation_band = dataset.GetRasterBand(1)

elevation_data = elevation_band.ReadAsArray()



# 获取地理信息

geotransform = dataset.GetGeoTransform()

origin_x = geotransform[0]

origin_y = geotransform[3]

pixel_width = geotransform[1]

pixel_height = geotransform[5]



# 定义网格参数

lon_min = 110.0

lon_max = 120.0

lat_min = 30.0

lat_max = 40.0

resolution = 0.1



# 计算网格尺寸

n_lon = int((lon_max - lon_min) / resolution) + 1

n_lat = int((lat_max - lat_min) / resolution) + 1



# 生成网格

lons = np.linspace(lon_min, lon_max, n_lon)

lats = np.linspace(lat_min, lat_max, n_lat)



# 将 DEM 数据与网格相结合

# 假设 DEM 数据的分辨率与网格分辨率相同

# 如果不同，需要进行插值处理

terrain_data = elevation_data



# 保存地形数据到文件

terrain_file = 'path/to/output/terrain.txt'

header = f'ncols {n_lon}\n'

header += f'nrows {n_lat}\n'

header += f'xllcorner {origin_x}\n'

header += f'llcorner {origin_y}\n'

header += f'cellsize {resolution}\n'

header += 'NODATA_value -9999\n'



with open(terrain_file, 'w') as f:

    f.write(header)

    np.savetxt(f, terrain_data, fmt='%1.2f')

5. 网格的可视化

可视化网格可以帮助我们更好地理解网格设置和地形数据。以下是一个示例代码，展示如何使用 Matplotlib 可视化水平网格和地形数据：


# 导入必要的库

import numpy as np

import matplotlib.pyplot as plt

from osgeo import gdal



# 读取 DEM 数据

dem_file = 'path/to/your/dem/file.dem'

dataset = gdal.Open(dem_file)

elevation_band = dataset.GetRasterBand(1)

elevation_data = elevation_band.ReadAsArray()



# 获取地理信息

geotransform = dataset.GetGeoTransform()

origin_x = geotransform[0]

origin_y = geotransform[3]

pixel_width = geotransform[1]

pixel_height = geotransform[5]



# 定义网格参数

lon_min = 110.0

lon_max = 120.0

lat_min = 30.0

lat_max = 40.0

resolution = 0.1



# 计算网格尺寸

n_lon = int((lon_max - lon_min) / resolution) + 1

n_lat = int((lat_max - lat_min) / resolution) + 1



# 生成网格

lons, lats = np.meshgrid(np.linspace(lon_min, lon_max, n_lon), np.linspace(lat_min, lat_max, n_lat))



# 可视化地形数据

plt.figure(figsize=(10, 10))

plt.contourf(lons, lats, elevation_data, cmap='terrain')

plt.colorbar(label='Elevation (m)')

plt.xlabel('Longitude (°)')

plt.ylabel('Latitude (°)')

plt.title('Terrain Elevation')

plt.show()

6. 网格设置的优化

合理的网格设置可以提高模型的精度和计算效率。以下是一些优化网格设置的建议：

6.1 选择合适的网格分辨率

高分辨率：适用于小范围、高精度的模拟。
低分辨率：适用于大范围、低精度的模拟。

6.2 选择合适的网格投影

经纬度投影：适用于大范围的模拟。
UTM 投影：适用于小范围的模拟，可以减少投影误差。

6.3 优化垂直网格设置

均匀分布：适用于简单的垂直结构模拟。
分层分布：适用于复杂的垂直结构模拟，如边界层和自由大气层。

7. 地形和网格设置的常见问题

在进行地形和网格设置时，可能会遇到一些常见问题。以下是一些解决方案：

7.1 地形数据分辨率不匹配

如果地形数据的分辨率与网格分辨率不匹配，可以使用插值方法进行处理。例如，使用 SciPy 库中的 griddata 函数进行插值：


# 导入必要的库

import numpy as np

import scipy.interpolate as interpolate

from osgeo import gdal



# 读取 DEM 数据

dem_file = 'path/to/your/dem/file.dem'

dataset = gdal.Open(dem_file)

elevation_band = dataset.GetRasterBand(1)

elevation_data = elevation_band.ReadAsArray()



# 获取地理信息

geotransform = dataset.GetGeoTransform()

origin_x = geotransform[0]

origin_y = geotransform[3]

pixel_width = geotransform[1]

pixel_height = geotransform[5]



# 定义目标网格参数

lon_min = 110.0

lon_max = 120.0

lat_min = 30.0

lat_max = 40.0

resolution = 0.1



# 计算目标网格尺寸

n_lon = int((lon_max - lon_min) / resolution) + 1

n_lat = int((lat_max - lat_min) / resolution) + 1



# 生成目标网格

target_lons, target_lats = np.meshgrid(np.linspace(lon_min, lon_max, n_lon), np.linspace(lat_min, lat_max, n_lat))



# 将 DEM 数据与目标网格相结合

# 使用插值方法

points = np.array([[i * pixel_width + origin_x, j * pixel_height + origin_y] for i in range(elevation_data.shape[1]) for j in range(elevation_data.shape[0])])

values = elevation_data.flatten()



# 插值

interpolated_elevation = interpolate.griddata(points, values, (target_lons, target_lats), method='linear')



# 保存插值后的地形数据到文件

terrain_file = 'path/to/output/interpolated_terrain.txt'

header = f'ncols {n_lon}\n'

header += f'nrows {n_lat}\n'

header += f'xllcorner {lon_min}\n'

header += f'llcorner {lat_min}\n'

header += f'cellsize {resolution}\n'

header += 'NODATA_value -9999\n'



with open(terrain_file, 'w') as f:

    f.write(header)

    np.savetxt(f, interpolated_elevation, fmt='%1.2f')

7.2 网格投影选择不当

如果选择的网格投影方式不当，可能会导致投影误差。建议在 GIS 软件中进行投影转换，确保数据的一致性。例如，使用 ArcGIS 的 Project 工具将数据从一种投影方式转换为另一种。

8. 实例应用

以下是一个完整的实例，展示如何从 USGS 下载 DEM 数据，处理并设置地形和网格参数，最后进行可视化。

8.1 下载 DEM 数据

假设我们已经从 USGS 下载了 DEM 数据文件 dem_file.dem。

8.2 处理地形数据


# 导入必要的库

import numpy as np

from osgeo import gdal

import scipy.interpolate as interpolate

import matplotlib.pyplot as plt



# 读取 DEM 数据

dem_file = 'path/to/your/dem/dem_file.dem'

dataset = gdal.Open(dem_file)

elevation_band = dataset.GetRasterBand(1)

elevation_data = elevation_band.ReadAsArray()



# 获取地理信息

geotransform = dataset.GetGeoTransform()

origin_x = geotransform[0]

origin_y = geotransform[3]

pixel_width = geotransform[1]

pixel_height = geotransform[5]



# 定义目标网格参数

lon_min = 110.0

lon_max = 120.0

lat_min = 30.0

lat_max = 40.0

resolution = 0.1



# 计算目标网格尺寸

n_lon = int((lon_max - lon_min) / resolution) + 1

n_lat = int((lat_max - lat_min) / resolution) + 1



# 生成目标网格

target_lons, target_lats = np.meshgrid(np.linspace(lon_min, lon_max, n_lon), np.linspace(lat_min, lat_max, n_lat))



# 将 DEM 数据与目标网格相结合

# 使用插值方法

points = np.array([[i * pixel# 地形和网格设置



## 地形数据的准备



在进行空气质量仿真时，地形数据是至关重要的输入之一。地形数据的准确性和分辨率直接影响到模型的模拟结果。CAMx 支持多种地形数据格式，常见的包括 DEM（Digital Elevation Model）和 TERR（Terrain Data）格式。地形数据通常需要从外部数据源获取，如 USGS（United States Geological Survey）或 GIS（Geographic Information System）软件。



### 1. 获取地形数据



#### 1.1 从 USGS 获取 DEM 数据



USGS 提供了大量的地形数据，可以通过其官方网站下载。以下是一个简单的步骤说明：



1. 访问 USGS 网站。

2. 选择你感兴趣的区域。

3. 下载所需的 DEM 数据文件，通常是 `.hgt` 或 `.dem` 格式。



#### 1.2 从 GIS 软件获取 TERR 数据



使用 GIS 软件如 ArcGIS，可以生成 TERR 格式的地形数据。以下是一个简单的步骤说明：



1. 打开 ArcGIS。

2. 添加你要使用的 DEM 数据。

3. 使用工具将 DEM 数据转换为 TERR 格式。例如，使用 `Raster to ASCII` 工具将 DEM 转换为 ASCII 格式，然后进一步转换为 TERR 格式。



### 2. 地形数据的处理



#### 2.1 使用 GDAL 工具处理地形数据



GDAL（Geospatial Data Abstraction Library）是一个强大的开源地理空间数据处理库，可以用于处理 DEM 数据。以下是一个示例代码，展示如何使用 GDAL 将 DEM 数据转换为 ASCII 格式：



```python

# 导入 GDAL 库

from osgeo import gdal

import numpy as np



# 读取 DEM 数据

dem_file = 'path/to/your/dem/file.dem'

dataset = gdal.Open(dem_file)



# 获取地形数据

elevation_band = dataset.GetRasterBand(1)

elevation_data = elevation_band.ReadAsArray()



# 获取地理信息

geotransform = dataset.GetGeoTransform()

origin_x = geotransform[0]

origin_y = geotransform[3]

pixel_width = geotransform[1]

pixel_height = geotransform[5]



# 写入 ASCII 文件

output_file = 'path/to/output/terrain.asc'

header = f'ncols {elevation_data.shape[1]}\n'

header += f'nrows {elevation_data.shape[0]}\n'

header += f'xllcorner {origin_x}\n'

header += f'llcorner {origin_y}\n'

header += f'cellsize {pixel_width}\n'

header += 'NODATA_value -9999\n'



with open(output_file, 'w') as f:

    f.write(header)

    np.savetxt(f, elevation_data, fmt='%1.2f')