轻磅！Python风场流线图与三种滤波方法

用户11172986

发布于 2024-06-20 18:02:13

890

发布于 2024-06-20 18:02:13

前言

感谢读者【摸鱼大佬】来信提供程序与文件九点平滑的工作原理是将风速数据中的每个值替换为该值及其八个相邻值的平均值。这具有平滑数据和消除任何高频噪声的效果。解释九点平滑器是如何工作的：创建一个新数组来存储平滑后的值。迭代风速数据数组。对于数组中的每个值，获取八个相邻值。计算九个值的平均值。在新数组中存储平均值。平滑过程完成后，新阵列将包含平滑后的风速数据。

下面测试九点平滑下的台风流线图与scipy的滤波结果进行对比核心函数：scipy.ndimage.gaussian_filter scipy.ndimage.median_filter

九点平滑

In [1]:

%time
"""
Created on Fri Jul 29 20:01:19 2022

@author: asus
"""

import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.mpl.ticker as cticker
import cartopy.io.shapereader as shpreader
import numpy as np
import matplotlib.patches as mpatches
plt.rcParams['font.sans-serif'] = ['Source Han Sans CN']  
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
#%% 九点滤波函数
def smooth9(array):
    a, b = array.shape
    
    new_array = array.copy()
    base_array = array.copy()
    
    for t in range(2):
        for i in range(1,a-1):
            for j in range(1,b-1):
                new_array[i,j]  = (base_array[i,j-1:j+2].sum() + base_array[i+1,j-1:j+2].sum() + base_array[i-1,j-1:j+2].sum())/9
        base_array = new_array.copy()
    return array - new_array

#%% 
for ti in ['2022-07-01-23','2022-07-02-01','2013-07-01-23','2013-07-02-01']:
    if ti[0:4] == '2013':
        ex = 0
        name = ti + '温比亚'
    else:
        ex = 1
        name = ti + '暹芭'
    # sst = xr.open_dataset('SST.nc')['sst'].sel(time=ti2,longitude=slice(90,150),latitude=slice(40,0))[0,ex,:,:] - 273.15
    data1 = xr.open_dataset('/home/mw/project/GUV850_ver2.nc').sel(longitude=slice(90,150),latitude=slice(40,0))
    gpm = data1['z'].sel(time=ti)[ex,:,:] / 9.8 / 10
    u = data1['u'].sel(time=ti)[ex,:,:]
    v = data1['v'].sel(time=ti)[ex,:,:]
    
    # 空间滤波
    u.values = smooth9(u.values)
    v.values = smooth9(v.values)
    lon, lat = np.meshgrid(data1.longitude, data1.latitude)
    
    #%% 
    fig = plt.figure(figsize=(13,8))#创建图窗
    proj = ccrs.PlateCarree(central_longitude=180)
    
    # 地图范围
    leftlon, rightlon, lowerlat, upperlat = [108.6,112.8,18,23]
        
    ax = fig.add_subplot(projection = proj)
    
    
    # ax.stock_img()
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'),edgecolor='grey')#海岸线
    ax.set_xticks(np.arange(leftlon+0.2,rightlon-0.2+0.5,0.5), crs=ccrs.PlateCarree())#经度标签
    ax.set_yticks(np.arange(lowerlat+0.2,upperlat-0.2+0.4,0.4), crs=ccrs.PlateCarree())#纬度标签
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    ax.set_extent([leftlon,rightlon,lowerlat,upperlat],crs=ccrs.PlateCarree())#设置图形
    
    
    # 读取地图文件
    china = shpreader.Reader('/home/mw/input/china1656/china_map/china_map/China_Province_2022.shp').geometries()
    ax.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='grey',linewidth=0.5,
                      zorder = 1)
    

    # 数据切割
    u = u.sel(longitude=slice(leftlon-1,rightlon+1),latitude=slice(upperlat+1,lowerlat-1))
    v = v.sel(longitude=slice(leftlon-1,rightlon+1),latitude=slice(upperlat+1,lowerlat-1))
    lon, lat = np.meshgrid(u.longitude, v.latitude)
    
    # 流线
    st = ax.streamplot(lon,lat,u.values,v.values,transform=ccrs.PlateCarree(),
                        density=3)
    # 标题
    ax.set_title(name,loc='left',fontdict={'size': 18})

CPU times: user 2 µs, sys: 0 ns, total: 2 µs
Wall time: 4.05 µs

/opt/conda/lib/python3.9/site-packages/xarray/backends/plugins.py:71: RuntimeWarning: Engine 'cfradial1' loading failed:
cannot import name 'ParasiteAxesAuxTrans' from 'mpl_toolkits.axisartist' (/opt/conda/lib/python3.9/site-packages/mpl_toolkits/axisartist/__init__.py)
  warnings.warn(f"Engine {name!r} loading failed:\n{ex}", RuntimeWarning)
/opt/conda/lib/python3.9/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

中值滤波

函数的参数介绍如下：

input：要进行中值滤波的输入图像或数组。 size：一个整数或整数元组，指定滤波器的大小。如果是整数，则表示每个维度上的滤波器大小相同。如果是元组，则可以为每个维度指定不同的滤波器大小。 footprint：指定滤波器的结构元素。它是一个布尔型数组，用于定义滤波器的形状。默认值为 None，表示使用一个与滤波器大小相同的全连接结构元素。 output：用于存储结果的数组。如果未提供，则会创建一个与输入数组相同类型和形状的新数组。 mode：边界模式。默认值为 'reflect'，表示对超出边界的像素进行镜像反射处理。其他可选值包括 'constant'（使用常数填充），'nearest'（使用最近的边界像素填充）和 'wrap'（循环填充）。 cval：当 mode 为 'constant' 时使用的常数值。 origin：滤波器应用的锚点位置。默认值为 -1，表示滤波器中心。 brute_force：是否使用暴力方法进行计算。默认值为 False，表示使用快速方法进行计算。该函数返回中值滤波后的图像或数组。

In [11]:

%time
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.ticker as cticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import shapely.geometry as sgeom
from concurrent.futures import ThreadPoolExecutor
from scipy.ndimage.filters import gaussian_filter, median_filter
# # 九点滤波函数
# def smooth9(array):
#     a, b = array.shape
#     new_array = array.copy()

#     for t in range(2):
#         new_array[1:-1, 1:-1] = (array[:-2, 1:-1] + array[1:-1, :-2] + array[1:-1, 1:-1] + array[1:-1, 2:] + array[2:, 1:-1]) / 5
#         array = new_array.copy()

#     return array

def plot_streamline(ti, ex):
    if ti[0:4] == '2013':
        ex = 0
        name = ti + '温比亚'
    else:
        ex = 1
        name = ti + '暹芭'
    
    data1 = xr.open_dataset('/home/mw/project/GUV850_ver2.nc').sel(longitude=slice(90,150),latitude=slice(40,0))
    gpm = data1['z'].sel(time=ti)[ex,:,:] / 9.8 / 10
    u = data1['u'].sel(time=ti)[ex,:,:]
    v = data1['v'].sel(time=ti)[ex,:,:]

    # 空间滤波
    u.values =median_filter(u.values , size=1)
    v.values =median_filter(v.values , size=1)
    # 数据切割
    leftlon, rightlon, lowerlat, upperlat = [108.6, 112.8, 18, 23]
    u = u.sel(longitude=slice(leftlon-1,rightlon+1),latitude=slice(upperlat+1,lowerlat-1))
    v = v.sel(longitude=slice(leftlon-1,rightlon+1),latitude=slice(upperlat+1,lowerlat-1))

    # 创建图窗
    fig = plt.figure(figsize=(13,8))
    proj = ccrs.PlateCarree(central_longitude=180)
    
    # 地图范围
    ax = fig.add_subplot(projection=proj)
    ax.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())

    # 添加地理特征
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'), edgecolor='grey')
    china = cfeature.ShapelyFeature(shpreader.Reader('/home/mw/input/china1656/china_map/china_map/China_Province_2022.shp').geometries(),
                                    ccrs.PlateCarree(), facecolor='none', edgecolor='grey', linewidth=0.5)
    ax.add_feature(china)

    # 绘制流线图
    lon, lat = np.meshgrid(u.longitude, v.latitude)
    st = ax.streamplot(lon, lat, u.values, v.values, transform=ccrs.PlateCarree(), density=3)

    # 设置标题
    ax.set_title(name, loc='left', fontdict={ 'size': 18})

    # 保存图像
    fig.savefig(ti + '.jpg', dpi=600)

# 并行计算绘制流线图
with ThreadPoolExecutor() as executor:
    ti_list = ['2022-07-01-23', '2022-07-02-01', '2013-07-01-23', '2013-07-02-01']
    ex_list = [1, 1, 0, 0]
    executor.map(plot_streamline, ti_list, ex_list)

CPU times: user 3 µs, sys: 1e+03 ns, total: 4 µs
Wall time: 5.25 µs

<ipython-input-11-0911adaaf8f6>:10: DeprecationWarning: Please use `gaussian_filter` from the `scipy.ndimage` namespace, the `scipy.ndimage.filters` namespace is deprecated.
  from scipy.ndimage.filters import gaussian_filter, median_filter
<ipython-input-11-0911adaaf8f6>:10: DeprecationWarning: Please use `median_filter` from the `scipy.ndimage` namespace, the `scipy.ndimage.filters` namespace is deprecated.
  from scipy.ndimage.filters import gaussian_filter, median_filter

高斯滤波

函数的参数介绍如下：

input：要进行高斯滤波的输入图像或数组。 sigma：高斯核的标准差，控制滤波的平滑程度。可以是一个标量，表示在每个维度上使用相同的标准差；也可以是一个具有与输入数组相同数量的元素的序列，分别指定每个维度上的标准差。 order：指定高斯滤波器的导数阶数。默认为 0，表示只进行平滑操作。如果设置为 1，则会计算一阶导数；如果设置为 2，则会计算二阶导数，以此类推。 output：用于存储结果的数组。如果未提供，则会创建一个与输入数组相同类型和形状的新数组。 mode：边界模式。默认值为 'reflect'，表示对超出边界的像素进行镜像反射处理。其他可选值包括 'constant'（使用常数填充），'nearest'（使用最近的边界像素填充）和 'wrap'（循环填充）。 cval：当 mode 为 'constant' 时使用的常数值。该函数返回经过高斯滤波后的图像或数组。

In [10]:

%time
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.ticker as cticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import shapely.geometry as sgeom
from concurrent.futures import ThreadPoolExecutor
from scipy.ndimage.filters import gaussian_filter, median_filter

def plot_streamline(ti, ex):
    if ti[0:4] == '2013':
        name = ti + '温比亚'
    else:
        name = ti + '暹芭'
    
    data1 = xr.open_dataset('/home/mw/project/GUV850_ver2.nc').sel(longitude=slice(90,150),latitude=slice(40,0))
    gpm = data1['z'].sel(time=ti)[ex,:,:] / 9.8 / 10
    u = data1['u'].sel(time=ti)[ex,:,:]
    v = data1['v'].sel(time=ti)[ex,:,:]

    # 空间滤波
    u.values =gaussian_filter(u.values , sigma=1)
    v.values =gaussian_filter(v.values , sigma=1)
    # 数据切割
    leftlon, rightlon, lowerlat, upperlat = [108.6, 112.8, 18, 23]
    u = u.sel(longitude=slice(leftlon-1,rightlon+1),latitude=slice(upperlat+1,lowerlat-1))
    v = v.sel(longitude=slice(leftlon-1,rightlon+1),latitude=slice(upperlat+1,lowerlat-1))

    # 创建图窗
    fig = plt.figure(figsize=(13,8))
    proj = ccrs.PlateCarree(central_longitude=180)
    
    # 地图范围
    ax = fig.add_subplot(projection=proj)
    ax.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())

    # 添加地理特征
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'), edgecolor='grey')
    china = cfeature.ShapelyFeature(shpreader.Reader('/home/mw/input/china1656/china_map/china_map/China_Province_2022.shp').geometries(),
                                    ccrs.PlateCarree(), facecolor='none', edgecolor='grey', linewidth=0.5)
    ax.add_feature(china)

    # 绘制流线图
    lon, lat = np.meshgrid(u.longitude, v.latitude)
    st = ax.streamplot(lon, lat, u.values, v.values, transform=ccrs.PlateCarree(), density=3)

    # 设置标题
    ax.set_title(name, loc='left', fontdict={ 'size': 18})

    # 保存图像
    fig.savefig(ti + '.jpg', dpi=600)

# 并行计算绘制流线图
with ThreadPoolExecutor() as executor:
    ti_list = ['2022-07-01-23', '2022-07-02-01', '2013-07-01-23', '2013-07-02-01']
    ex_list = [1, 1, 0, 0]
    executor.map(plot_streamline, ti_list, ex_list)

CPU times: user 4 µs, sys: 1e+03 ns, total: 5 µs
Wall time: 6.91 µs

<ipython-input-10-0f2f7fb0d73e>:10: DeprecationWarning: Please use `gaussian_filter` from the `scipy.ndimage` namespace, the `scipy.ndimage.filters` namespace is deprecated.
  from scipy.ndimage.filters import gaussian_filter, median_filter
<ipython-input-10-0f2f7fb0d73e>:10: DeprecationWarning: Please use `median_filter` from the `scipy.ndimage` namespace, the `scipy.ndimage.filters` namespace is deprecated.
  from scipy.ndimage.filters import gaussian_filter, median_filter