气象编程 | 数据筛选1-站点数据筛选

站点数据筛选

站点数据是常规观测中一种常见的气象数据保存样式，一般是基于站点实况的，对于这种数据的筛选与白化方法是多种多样的，这一章节讨论一些常见的站点数据筛选与白化方法。

前期准备

导入一些常见的站点数据处理与可视化方法库包，基本设置等，这里我们选用了一份电子表格样式的站点数据，数据被人为添加随机抖动加密：

import pandas as pd
import numpy as np
import xarray as xr
import rioxarray as rioxr
import geopandas as gpd
import cartopy.crs as ccrs
import cartopy.io.shapereader as cis
import cartopy.mpl.ticker as cmt
import cartopy.mpl.patch as cmp
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.colors as mcolors
import matplotlib.path as mpath
plt.rcParams['font.sans-serif']=['FangSong']
plt.rcParams['axes.unicode_minus']=False
file=r'C:\Desktop\天气2024年05月14日14时.xls'

读取电子表格文件，并添加随机抖动加密，读取地图文件。

df=pd.read_excel(file,skiprows=[0],na_values=['-','----','////',999999])
df['经度']=df['经度']+np.random.rand(len(df['经度']))*0.1
df['纬度']=df['纬度']+np.random.rand(len(df['纬度']))*0.1
gdf=gpd.read_file(r'F:\shp文件合集\恩施土家族苗族自治州_行政边界\恩施土家族苗族自治州_行政边界.shp',encoding='gbk')

电子表格数据格式：

通过行政区划筛选数据

如果数据格式有行政区划列，则可以通过行政区划列方便的提取所需数据，对于单个区划可以使用求等符号==，对于多个行政区划，可以使用isin方法。

df[df['县']=='利川市']

df[df['县'].isin(['利川市','宣恩县'])]

通过geopandas筛选数据

大多数时候的站点数据，都没有明确的行政区划列，但一般会携带经纬度列，可以通过经纬度数据来实现数据筛选。
首先将数据文件中的经纬度列转化为geometry。geopandas提供了快速命令from_xy。要不要设置crs看实际需求。epsg:4326就是cartopy的PlateCarree。

gpd.GeoSeries.from_xy(df['经度'],df['纬度']).set_crs('epsg:4326')

随后将geometry添加到原本的数据中，将数据样式转化为geopandas.GeoDataFrame。此时在数据最后新增一列geometry。

df_geo=gpd.GeoDataFrame(df,geometry=gpd.GeoSeries.from_xy(df['经度'],df['纬度']).set_crs('epsg:4326'))

读取我们需要的目标区县几何数据，这里还是恩施全州8个区县的信息： 仅提取利川市：

gdf[gdf['Name']=='利川市']

随后使用clip命令裁剪站点数据：

df_geo.clip(gdf[gdf['Name'




    
]=='利川市'])

先绘图再白化

如果仅需要绘图，可以通过先绘图再裁剪白化的方式实现数据按区划绘制。

points=ax4.scatter(df['经度'],df['纬度'],s=5,zorder=6)
points.set_clip_path(mpath.Path.make_compound_path(*cmp.geos_to_path(gdf[gdf['Name']=='利川市'].geometry.iloc[0])),
                     transform=ax4.transData)

各种方法绘图效果

fig=plt.figure(figsize=(6,6),dpi=700)
ax1=plt.subplot(221,projection=ccrs.PlateCarree())
ax2=plt.subplot(222,projection=ccrs.PlateCarree())
ax3=plt.subplot(223,projection=ccrs.PlateCarree())
ax4=plt.subplot(224,projection=ccrs.PlateCarree())
def draw_map(ax):
    ax.add_geometries(cis.Reader(r'E:\shp\恩施土家族苗族自治州_行政边界\恩施土家族苗族自治州_行政边界.shp').geometries(),
                      crs=ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=0.5)
    ax.add_geometries(gdf[gdf['Name']=='利川市'].geometry,
                      crs=ccrs.PlateCarree(),facecolor='none',edgecolor='tab:red',linewidth=1)
    ax.set_extent([108,111,28.75,31.75])
    ax.set_xticks(np.arange(108,111.1,1))
    ax.set_yticks(np.arange(28.75,31.76,1))
    ax.tick_params(labelsize=7)
    ax.xaxis.set_major_formatter(cmt.LongitudeFormatter())
    ax.yaxis.set_major_formatter(cmt.LatitudeFormatter())
[draw_map(ax) for ax in [ax1,ax2,ax3,ax4]]
[ax.set_xticklabels([]) for ax in [ax1,ax2]]
[ax.set_yticklabels([]) for ax in [ax2,ax4]]
ax1.set_title('全部站点分布',fontsize=12)
ax2.set_title('区划筛选站点分布',fontsize=12)
ax3.set_title('GeoPandas筛选站点分布',fontsize=12)
ax4.set_title('全部站点绘制后白化',fontsize=12)
ax1.scatter(df['经度'],df['纬度'],s=5,zorder=6)
ax2.scatter(df[df['县']=='利川市']['经度'],df[df['县']=='利川市']['纬度'],s=5,zorder=6)
ax3.scatter(df_geo.clip(gdf[gdf['Name']=='利川市'])['经度'],df_geo.clip(gdf[gdf['Name']=='利川市'])['纬度'],s=5,zorder=6)
points=ax4.scatter(df['经度'],df['纬度'],s=5,zorder=6)
points.set_clip_path(mpath.Path.make_compound_path(*cmp.geos_to_path(gdf[gdf['Name']=='利川市'].geometry.iloc[0])),
                     transform=ax4.transData)

请注意，绘图后再set_clip_path是真正的裁剪操作，会将图形按照给出的路径裁剪： 下一章节为栅格数据的筛选。