Python气象数据处理与可视化

import cdsapi
import os
from subprocess import call


def idmDownloader(task_url, folder_path, file_name):
    """
    IDM下载器
    :param task_url: 下载任务地址
    :param folder_path: 存放文件夹
    :param file_name: 文件名
    :return:
    """
    # IDM安装目录
    idm_engine = "C:\\Program Files (x86)\\Internet Download Manager\\IDMan.exe"
    # 将任务添加至队列
    call([idm_engine, '/d', task_url, '/p', folder_path, '/f', file_name, '/a'])
    # 开始任务队列
    call([idm_engine, '/s'])


c = cdsapi.Client()

for date_str in extreme_dates:
    year, month, day = date_str.split('-')
    year_path = os.path.join(path, year)  # 创建每年的文件夹
    os.makedirs(year_path, exist_ok=True)  # 如果文件夹不存在则创建

    filename = f'{variable}.0.25deg.{year}.{month}.{day}.nc'
    output_file_path = os.path.join(year_path, filename)

    # 检查文件是否已存在
    if os.path.exists(output_file_path):
        print(f'文件 {output_file_path} 已存在，跳过下载。')
        continue

    print(f'=======  日期:  {date_str}  =======')

    r = c.retrieve(
        'derived-era5-pressure-levels-daily-statistics',
        {
            'product_type': 'reanalysis',
            'variable': ["geopotential",
                         "specific_humidity",
                         "u_component_of_wind",
                         "v_component_of_wind",
                         "vertical_velocity"],
            'year': [year],
            'month': [month],
            'day': [day],
            "pressure_level": ["500"],
            "daily_statistic": "daily_mean",
            "time_zone": "utc+08:00",
            "frequency": "1_hourly",
            "area": [80, 30, 5, 150]
        },
    )

    url = r.location  # 获取文件下载地址
    idmDownloader(url, year_path, filename)  # 添加进IDM中下载

https://disc.gsfc.nasa.gov/

import xarray as xr

# 打开文件
ds = xr.open_dataset('example.nc')

# 查看数据集的基本信息
print(ds)

# 读取某个变量的数据
temperature = ds['temperature'].values
print("Temperature Data:", temperature)

# 关闭数据集
ds.close()

import netCDF4 as nc

# 打开文件
dataset = nc.Dataset('example.nc', 'r')

# 查看文件中的维度和变量
print("Dimensions:", dataset.dimensions.keys())
print("Variables:", dataset.variables.keys())

# 读取某个变量的数据
temperature = dataset.variables['temperature'][:]
print("Temperature Data:", temperature)

# 关闭文件
dataset.close()

import netCDF4 as nc

# 打开NetCDF文件
dataset = nc.Dataset('example.nc', 'r')

# 查看文件中的变量
print("Variables:", dataset.variables.keys())

# 索引并读取特定变量的数据
temperature = dataset.variables['temperature'][:]  # 读取所有数据
temperature_at_time_0 = dataset.variables['temperature'][0, :, :]  # 索引时间维度的第一个时间点
temperature_at_lat_0_lon_0 = dataset.variables['temperature'][:, 0, 0]  # 索引特定经纬度点

print("Temperature Data:", temperature)
print("Temperature at Time 0:", temperature_at_time_0)
print("Temperature at Lat 0, Lon 0:", temperature_at_lat_0_lon_0)

# 关闭文件
dataset.close()

import xarray as xr

# 打开NetCDF文件
ds = xr.open_dataset('example.nc')

# 查看数据集信息
print(ds)

# 索引特定变量的数据
temperature = ds['temperature']  # 获取温度变量
temperature_at_time_0 = temperature.isel(time=0)  # 索引时间维度的第一个时间点
temperature_at_lat_0_lon_0 = temperature.isel(lat=0, lon=0)  # 索引特定经纬度点

print("Temperature Data:", temperature)
print("Temperature at Time 0:", temperature_at_time_0)
print("Temperature at Lat 0, Lon 0:", temperature_at_lat_0_lon_0)

# 关闭数据集
ds.close()

slice(start, stop, step)

import rioxarray
import geopandas
from shapely. geometry import mapping
#设置数据集的空间维度（也称为地理或空间坐标）并编写坐标参考系统（crs）。 EPSG4326 坐标系，也称为 WGS 84（世界大地测量系统 1984），是地理数据常用的坐标系
data.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True)    # x_dim具体是什么，看自己的nc文件变量名
data.rio.write_crs("epsg:4326", inplace=True)
# 将非空间变量设置为坐标，避免干扰裁剪操作
data = data.set_coords(["time_bnds"])  # 如果nc文件里有lon_bnds啥的也要加上

hunan = geopandas.read_file('C:/Users/lenovo/Desktop/湖南无子区域/湖南省行政边界.shp', crs="epsg:4326")  # 该shp只含湖南省界
clipped = data.rio.clip(hunan.geometry.apply(mapping), hunan.crs, drop=True)

print(clipped)

#coding=gbk
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata

# 1. 打开数据
file_path = r'C:/path/to/your/data.nc'  # 替换为你的NC文件路径
ds = xr.open_dataset(file_path)

# 打印数据集信息
print(ds)

# 2. 索引数据
# 假设数据集中有一个变量名为'temperature'，维度为(time, lat, lon)
temperature = ds['temperature']  # 提取温度变量
print(temperature)

# 3. 裁剪数据
# 假设需要裁剪到特定区域（如纬度20°-40°，经度100°-120°）
lat_range = slice(20, 40)  # 纬度范围
lon_range = slice(100, 120)  # 经度范围
temperature_cropped = temperature.sel(lat=lat_range, lon=lon_range)
print(temperature_cropped)

# 4. 白化（掩膜处理）
# 假设需要掩膜掉海洋区域（海洋值为NaN）
# 这里假设海洋区域的温度值为NaN
temperature_masked = temperature_cropped.where(~np.isnan(temperature_cropped))
print(temperature_masked)

# 5. 插值
# 假设需要将数据插值到更精细的网格（如0.1°×0.1°）
# 原始网格
lon_original = temperature_masked.lon.values
lat_original = temperature_masked.lat.values

# 新网格
new_lon = np.linspace(lon_original.min(), lon_original.max(), int((lon_original.max() - lon_original.min()) / 0.1))
new_lat = np.linspace(lat_original.min(), lat_original.max(), int((lat_original.max() - lat_original.min()) / 0.1))
new_lon, new_lat = np.meshgrid(new_lon, new_lat)

# 插值（使用线性插值）
temperature_interpolated = griddata(
    (lon_original.flatten(), lat_original.flatten()),  # 原始网格点
    temperature_masked.values.flatten(),  # 原始数据
    (new_lon, new_lat),  # 新网格点
    method='linear'  # 插值方法：linear, nearest, cubic
)

# 6. 可视化结果
plt.figure(figsize=(10, 6))
plt.contourf(new_lon, new_lat, temperature_interpolated, levels=20, cmap='coolwarm')
plt.colorbar(label='Temperature (°C)')
plt.title('Interpolated Temperature')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.show()

# 7. 保存处理后的数据（可选）
output_path = r'C:/path/to/your/processed_data.nc'  # 替换为输出路径
temperature_interpolated_ds = xr.Dataset(
    {
        'temperature': (['lat', 'lon'], temperature_interpolated)
    },
    coords={
        'lon': new_lon[0, :],  # 新经度坐标
        'lat': new_lat[:, 0]  # 新纬度坐标
    }
)
temperature_interpolated_ds.to_netcdf(output_path)
print(f"处理后的数据已保存至：{output_path}")

fig, axes = plt.subplots(2, 2, figsize=(10, 8))
axes[0, 0].plot(x, np.sin(x), label="sin(x)")
axes[0, 1].plot(x, np.cos(x), label="cos(x)")
axes[1, 0].plot(x, np.tan(x), label="tan(x)")
axes[1, 1].plot(x, np.exp(-x), label="exp(-x)")
for ax in axes.flatten():
    ax.legend()
    ax.grid(True)
plt.tight_layout()
plt.show()

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

# 创建地图
fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={"projection": ccrs.PlateCarree()})

# 添加地图元素
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.add_feature(cfeature.LAND, facecolor="lightgray")
ax.add_feature(cfeature.OCEAN, facecolor="lightblue")

# 设置地图范围
ax.set_extent([70, 140, 15, 55], crs=ccrs.PlateCarree())  # 经纬度范围

# 添加标题
ax.set_title("China Map", fontsize=16)

plt.show()

import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

# 加载NetCDF数据
ds = xr.open_dataset("example.nc")
temperature = ds['temperature'].isel(time=0)

# 创建地图
fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={"projection": ccrs.PlateCarree()})

# 添加地图元素
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.add_feature(cfeature.LAND, facecolor="lightgray")

# 绘制数据
temperature.plot(ax=ax, transform=ccrs.PlateCarree(), cmap="coolwarm", add_colorbar=True)

# 设置地图范围
ax.set_extent([70, 140, 15, 55], crs=ccrs.PlateCarree())

# 添加标题
ax.set_title("Temperature Distribution", fontsize=16)

plt.show()