【WRF地图投影】详细解释WRF投影原理及操作步骤

【WRF地图投影】详细解释WRF投影原理及操作步骤

• • 坐标系统概述

• Map projections in WRF: it's easy to get it wrong

• • 如何绘制WRF的输出结果图？
  • WRF投影原理

• 如何将其应用于我们的 WRF 文件？

• • 情形1：使用salem包绘制WRF结果图
  • 情形2：自主确定投影并绘制图形
  • 基于cartopy包绘制

• 另：报错总结

• 参考

在学习WRF的过程中，偶尔看到一篇文章-Map projections in WRF: it’s easy to get it wrong。

由于WRF使用的是Lambert投影（公里坐标系），而我们一般使用栅格数据时都是用的 经纬度坐标系 ，有时也许会出现疑惑甚至是错误。因此本博客主要学习上述文章，并结合本人现有的对WRF模拟的过程进行一个解释和总结。如有任何错误或者建议，直接评论此博客，或者私信我都可。

坐标系统概述

坐标系统分为两类：
一、地理 坐标系统（Geographic Coordinate systems, GCS）：基于三维球体地球模型和经纬度坐标
二、投影 坐标系统（Projected Coordinate systems, PCS）：基于一个二维坐标平面，从球体地球的数学转换位置得到
两种类型的坐标系统都使用由球体的尺寸和球体与地球的匹配度（方向和原点）组成的基准。

常用的GCS坐标系统：WGS_1984（以 6378137.0m 为长半轴的椭球）
projection=regular_ll 表示规则的经纬度投影（Regular Latitude-Longitude Projection）

常用的PCS坐标系统：
projection = albers_nad83，在 NAD83 地理基准坐标系统 (GCS) 上使用了 Albers Equal-Area 投影（阿尔伯斯等面积投影），适用于北美
projection = ‘albers_cgcs2000’，基准坐标系通常使用 CGCS2000（中国大地坐标系 2000） 或 WGS84，适用于中国

说明：gis4wrf 仅支持基于球体（sphere）的投影基准面 ，且球体半径必须为 6370000m（WRF 默认的地球半径）

虽然同化到有限区域大气模式（LAM）中的许多地球物理数据集 是基于球体（即椭圆）地球模型，但大多数有限区域大气模式（LAM）假设地球表面是一个完美的球体，以简化计算。在中纬度地区，一个点的位置在球体（spheroidal Earth models）地球模型和球形（spherical Earth models）地球模型之间可能相差20公里以上（下图）。

Map projections in WRF: it’s easy to get it wrong

有关WRF Projection的总结如下：

• 您需要了解的有关 WRF 投影参数的所有信息都可以在 WRF netcdf 文件输出中找到
• 借助 WRF-Python 或 salem 包，您可以非常轻松地使用 python 绘制 WRF 输出
• 如果您需要自己解析投影参数：请不要忘记将 semimajor_axis 和 semiminor_axis 投影参数调整为 6370000（球体）

如何绘制WRF的输出结果图？

事实上，WRF 的地图投影非常特殊，以至于 NCAR 的一些科学家就此写了一篇论文。

论文-J2013-Overlapping Interests: The Impact of Geographic Coordinate
Assumptions on Limited-Area Atmospheric Model Simulations
下载链接-Monthly Weather Review

在高分辨率（10公里）应用中越来越多地使用有限区域模型（limited-area models, LAMs），为此，输入陆地（地形和土地利用数据）和气象数据集的一致映射 对于精确模拟至关重要。
大多数输入数据集的地理坐标系都是基于球形（即椭圆）地球模型，而LAMs通常假设地球是完美的球形。在预处理过程中，当输入数据被重新映射到LAM域时，这种区别经常被忽略，导致中纬度地区的地理位置差异可能超过20公里。
高分辨率LAM模拟对地球模式假设的敏感性强调了用户在预处理过程中确保地面和气象输入数据 一致映射的重要性（即，在重新映射到LAM域之前，数据集共享一个共同的地理坐标系 ）。同时，建模社区应该更新预处理系统，以确保为所有全局和有限区域的仿真域正确映射输入数据。

他们写道，地理定位错误可能导致“[…] 中纬度地区可能超过 20 公里的差异。[…] 同时，建模社区应该更新预处理系统，以确保输入数据正确映射到所有全球和有限区域模拟域。”

论文表明，这些差异可能会产生巨大影响。然而，这篇论文乍一看可能会产生误导：WRF 用户可以将其解读为计算标准 WGS84 椭球 和模型使用的球形地球之间的 基准偏移（Geodetic datum） 的动机。
事实上，你不应该这样做：WRF 默认地形和土地覆盖数据都在标准 WGS84 基准中提供。
这篇论文（以及另一篇令人惊讶的相似论文-Analysis of errors introduced by geographic coordinate systems on weather numeric prediction modeling）实际上认为，输入数据在被 WRF 使用之前也应该在球面坐标系中进行变换。我认为这会让事情变得更加复杂，所以我们最好暂时把这个问题放在一边。

如果你想现在停止阅读，只需记住这一点：
使用 WRF 球形地球半径参数来定义 WRF 地图投影，并将以这种方式获得的经度/纬度坐标视为在 WGS84 中有效。

WRF投影原理

该图说明了三点：
1、两个地图投影之间的经典转换需要三个步骤：
(1) 将东距/北距转换为经度/纬度
(2) 基准偏移（datum shift）
(3) 将经度/纬度转换回东距/北距

2、您必须使用 WRF 球面的半径 来定义 WRF 地图投影（在此示例中为兰伯特共形圆锥投影，但对于墨卡托或极地投影也是如此）

3、从球面 WRF 东距/北距转换中获得的经度/纬度坐标需要解释为标准 WGS84 基准中的经度/纬度，而不是球面经度/纬度。因此不需要基准偏移！

如何将其应用于我们的 WRF 文件？

情形1：使用salem包绘制WRF结果图

使用以下代码绘制geo_em.d01.nc文件，代码如下：

    import salem
    import matplotlib.pyplot as plt
    
    # 导入 NetCDF 数据
    file_path = 'D:/6 Python Codes/20241024 WRF Projection/geo_em.d01.nc'
    ds = salem.open_wrf_dataset(file_path)
    hgt_map = ds.HGT_M.where(ds.LANDMASK).salem.quick_map(cmap='topo');
    
    # 显示图像
    plt.show()
    
    # 关闭数据集
    ds.close()

绘制图形如下：

情形2：自主确定投影并绘制图形

现在，尽管 salem 可以为您完成所有工作，但让我们尝试在没有外部帮助的情况下正确进行投影。投影参数作为属性存储在 NetCDF 文件中：

    import pyproj
    wrf_proj = pyproj.Proj(proj='lcc', # projection type: Lambert Conformal Conic
                       lat_1=ds.TRUELAT1, lat_2=ds.TRUELAT2, # Cone intersects with the sphere
                       lat_0=ds.MOAD_CEN_LAT, lon_0=ds.STAND_LON, # Center point
                       a=6370000, b=6370000) # This is it! The Earth is a perfect sphere

现在，我们如何知道我们做对了？最好的测试是看看我们是否能够重新计算 WRF 在原始文件中提供的经度和纬度 ：

    # Easting and Northings of the domains center point
    wgs_proj = pyproj.Proj(proj='latlong', datum='WGS84')
    e, n = pyproj.transform(wgs_proj, wrf_proj, ds.CEN_LON, ds.CEN_LAT)
    # Grid parameters
    dx, dy = ds.DX, ds.DY
    nx, ny = ds.dims['west_east'], ds.dims['south_north']
    # Down left corner of the domain
    x0 = -(nx-1) / 2. * dx + e
    y0 = -(ny-1) / 2. * dy + n
    # 2d grid
    xx, yy = np.meshgrid(np.arange(nx) * dx + x0, np.arange(ny) * dy + y0)

xx 和 yy 是每个 WRF 网格点的二维坐标，以米为单位（WRF 地图投影中的东向和北向）。
第一个 WRF 域不需要计算域的中心点（它始终位于其自身投影的中心），但子域需要计算（子域可以放置在地图中的任何位置）。让我们将它们转换为经度和纬度，看看它们与原始 WRF 的比较结果如何：

    # Transform and plot
    our_lons, our_lats = pyproj.transform(wrf_proj, wgs_proj, xx, yy)
    ds['DIFF'] = np.sqrt((our_lons - ds.XLONG_M)**2 + (our_lats - ds.XLAT_M)**2)
    ds.salem.quick_map('DIFF', cmap='Reds');

我们的经度和纬度与 WRF 的逗号后最多 5 位数字一致。现在，如果我们忘记指定我们位于球体上，而是使用默认的 wgs84 椭球体，会发生什么情况？

    bad_proj = pyproj.Proj(proj='lcc', # projection type: Lambert Conformal Conic
                       lat_1=ds.TRUELAT1, lat_2=ds.TRUELAT2, # Cone intersects with the sphere
                       lat_0=ds.MOAD_CEN_LAT, lon_0=ds.STAND_LON, # Center point
                       ) # The Earth is now an ellipsoid
    bad_lons, bad_lats = pyproj.transform(bad_proj, wgs_proj, xx, yy)
    ds['DIFF2'] = np.sqrt((bad_lons - ds.XLONG_M)**2 + (bad_lats - ds.XLAT_M)**2)
    print('The max diff is: {}'.format(ds['DIFF2'].max().values))
    The max diff is: 0.114...

误差超过 0.1°！以米为单位计算误差更能说明问题。我们只需反转该过程：

    bad_xx, bad_yy = pyproj.transform(wgs_proj, bad_proj, ds.XLONG_M.values, ds.XLAT_M.values)
    ds['DIFF_M'] = np.sqrt((bad_xx - xx)**2 + (bad_yy - yy)**2) + ds.XLONG_M*0 # trick
    ds.salem.quick_map('DIFF_M', cmap='Reds');

地理定位中的误差在中心接近于零，而在域的上角上升到 7 公里。

好的……但我们为什么要关心呢？
如果您仅使用 WRF 提供的经度和纬度 进行分析，则不必关心这一点。它们在 WGS84 基准中有效，因此可用于查找距离气象站最近的网格点。

对于所有其他应用程序，您应该关心。

基于cartopy包绘制

    import cartopy
    import cartopy.crs as ccrs
    
    # Define the projection
    globe = ccrs.Globe(ellipse='sphere', semimajor_axis=6370000, semiminor_axis=6370000)
    lcc = ccrs.LambertConformal(globe=globe, # important!
                            central_longitude=ds.STAND_LON, central_latitude=ds.MOAD_CEN_LAT,
                            standard_parallels=(ds.TRUELAT1, ds.TRUELAT2),
                            )
    ax = plt.axes(projection=lcc)
    z.plot(ax=ax, transform=lcc, cmap='terrain');
    ax.coastlines()
    ax.add_feature(cartopy.feature.BORDERS, linestyle='-');
    ax.set_extent([xx.min(), xx.max(), yy.min(), yy.max()], crs=lcc)

另：报错总结

在使用转化器（transformer）时，出现以下错误：

    pyproj.exceptions.CRSError: Input is not a CRS: CONVERSION["PROJ-based coordinate operation",METHOD["PROJ-based operation method: +proj=longlat +datum=WGS84 +no_defs"]]

问题分析：当前的geopandas的版本太高了，函数失效引发的错误，降低该版本即可
解决办法：conda install geopandas==0.6.1

查看当前环境中版本geopandas，如下：

    conda activate myenv3.9
    conda list geopandas

(担心改变版本后，其他库包用不了，为了解决后患，先复制此环境 ，然后在新环境中尝试此解决方案)

    conda create -n myenv3.9.1 --clone myenv3.9
    conda install geopandas==0.6.1

参考

1、论文-J2013-Overlapping Interests: The Impact of Geographic Coordinate Assumptions on Limited-Area Atmospheric Model Simulations
下载链接-Monthly Weather Review
2、论文-J2017-Analysis of errors introduced by geographic coordinate systems on weather numeric prediction modeling
下载链接-EGU-Analysis of errors introduced by geographic coordinate systems on weather numeric prediction modeling
3、博客-pyproj.exceptions.CRSError: Invalid projection: +init=epsg:4326 +type=crs: (Internal Proj Error: pro