数据处理 | 投影转换案例：将非经纬度坐标数据转换至WGS84坐标系

用户11172986

发布于 2025-01-13 12:29:25

14500

代码可运行

文章被收录于专栏：气python风雨气python风雨

运行总次数：0

代码可运行

投影转换案例：将非经纬度坐标数据转换至WGS84坐标系

摘要

作者：崔忠强

编辑：气ython风雨

很多文件带着一些投影信息，导致经纬度和实际对应不上，这里提供一个做投影转换的案例。本文展示了如何读取这样的文件，并将其转换为常用的WGS84坐标系下的经纬度，以便进行正确的可视化和分析。

数据读取与初步检查

使用xarray库读取一个包含中国相对湿度月平均值的数据集文件。发现文件内部没有直接提供经纬度信息，只有x和y坐标以及时间维度。

import xarray as xr

ds = xr.open_dataset(r'path_to_file\HiTMC_Monthly_China_RH_201901_201912.nc')
print(ds)

获取投影信息并进行转换

`Projected Coordinate System: Albers_Conic_Equal_Area 
Projection: Albers 
false_easting: 4000000.00000000 
false_northing: 0.00000000 
central_meridian: 105.00000000 
standard_parallel_1: 25.00000000 
standard_parallel_2: 47.00000000 
latitude_of_origin: 0.00000000 
Linear Unit: Meter 
Geographic Coordinate System: GCS_WGS_1984 
Datum: D_WGS_1984 
Angular Unit: Degree `

根据附带的数据说明文档得知，该数据集使用的是阿尔伯斯等面积圆锥投影。因此，我们需要利用cartopy库来定义这个投影，并将其转换到WGS84坐标系下。

import cartopy.crs as ccrs
import numpy as np

# 构建网格
x, y = ds.x, ds.y
x_grid, y_grid = np.meshgrid(x.values, y.values)

# 定义源投影（阿尔伯斯）
albers_proj = ccrs.AlbersEqualArea(
    central_longitude=105.0,
    false_easting=4000000.0,
    false_northing=0.0,
    standard_parallels=(25.0, 47.0)
)

# 目标投影（WGS84）
wgs84 = ccrs.PlateCarree()

# 将格点化的x和y转换为一个个点坐标。
points = wgs84.transform_points(albers_proj, x_grid, y_grid)
lon, lat = points[:,:,0], points[:,:,1]
#这样我们就得到了所对应的经度和纬度了！
lon,lat

数据预处理

由于数据集中RH变量的单位是0.01%，需要进行适当的转换，并且替换掉异常值以确保数据的有效性。

variable_data = ds.RH.values / 100  # 单位转换为百分比
variable_data = variable_data[0,:,:]  # 选择1月份的数据
variable_data[variable_data > 100] = 100  # 剔除高异常值
variable_data[variable_data < 0] = 0  # 剔除低异常值或默认空值
variable_data

数据可视化

完成坐标转换和数据预处理后，可以使用matplotlib结合cartopy来进行数据的可视化展示。

import matplotlib.pyplot as plt
import cartopy.feature as cfeature

fig, ax = plt.subplots(subplot_kw={'projection': wgs84})
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.add_feature(cfeature.RIVERS)
ax.add_feature(cfeature.STATES)

contour = ax.contourf(lon, lat, variable_data, transform=wgs84, cmap='rainbow', vmin=0, vmax=100)
cbar = plt.colorbar(contour, orientation='horizontal', pad=0.05)
cbar.set_label('Relative Humidity (%)')  # 设置颜色条标签

ax.set_title('January 2019 Relative Humidity in WGS84 Coordinates')
plt.show()