序：
这篇博客旨在对赛题数据做一些初步的探索，包括查看数据中的缺失值、异常值等，以及通过可视化来观察各个特征的分布情况，为之后进行特征工程提供一些思路。

1. 查看数据整体情况

对于一份数据集，首先要对它的整体情况做一些基本的了解。

导入数据库

首先导入必要的库：

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns

from tqdm import tqdm
import multiprocessing as mp
import os
import pickle
import random

读取数据集

这份数据集是以压缩文件的形式给出的，其中包含每条渔船的csv格式的数据文件，我们这里定义读取数据文件的函数：

def read_train_file(filename=None):
    # 替换数据存放的路径
    Path = "../input/wisdomoceans/hy_round1_train_20200102/hy_round1_train_20200102/"
    return pd.read_csv(Path + filename,encoding="utf-8")

def read_test_file(filename=None):
    # 替换数据存放的路径
    Path = "../input/wisdomoceans/hy_round1_testA_20200102/hy_round1_testA_20200102/"
    return pd.read_csv(Path + filename,encoding="utf-8")

定义一个读取和保存数据的类：

# 定义加载和存储数据的类
class Load_Save_Data():
    def __init__(self, file_name=None):
        self.filename = file_name
        
    def load_data(self, Path=None):
        if Path is None:
            assert self.filename is not None, 'Invalid Path...'
        else:
            self.filename = Path
        with open(self.filename, 'wb') as f:
            data = pickle.load(f)
        return data
    
    def save_data(self, data, path):
        if path is None:
            assert self.filename is not None, 'Invalid path...'
        else:
            self.filename = path
        with open(self.filename, 'wb') as f:
            pickle.dump(data, f)

利用上述构建的函数和类，定义一个通过多线程的方式来读取数据的函数，加快读取速度：

def read_data(Path, Kind=''):
    filenames = os.listdir(Path)
    print('\n@Read Data From' + Path + '.......................')
    # 多线程处理
    with mp.Pool(processes=mp.cpu_count()) as pool:
        data_total = list(tqdm(pool.map(read_train_file if Kind=='train' else
                                        read_test_file, filenames), total=len(filenames)))
    print('\n@End Read Total Data .............................')
    load_save = Load_Save_Data()
    if Kind == 'train':
        load_save.save_data(data_total, './total_data.pkl')
    return data_total

现在可以用上面的函数来读取训练集和测试集：

# 读取训练数据
train_path = '../input/wisdomoceans/hy_round1_train_20200102/hy_round1_train_20200102/'
data_train = read_data(train_path, 'train')
data_train = pd.concat(data_train)

# 读取测试数据
test_path = '../input/wisdomoceans/hy_round1_testA_20200102/hy_round1_testA_20200102/'
data_test = read_data(test_path, 'test')
data_test = pd.concat(data_test)

查看数据集

先来看看训练集的基本情况：

print('The shape of train data:', data_train.shape)
data_train.head()

The shape of train data: (2699638, 7)

训练集的特征并不复杂，每条数据包括渔船ID、x轴坐标、y轴坐标、速度、方向、上报时间六个特征，要预测的目标值就是作业类型（type）。
查看训练集的各个特征的基本统计值：

data_train.describe([0.01, 0.025, 0.05, 0.5, 0.75, 0.9, 0.99])

可以看到速度的最大值与99%分位数相差很大，说明速度特征存在异常值，可以通过3-sigma原理判断异常值并进行处理。
最后查看训练集的缺失值：

# 查看训练集缺失值
print(f'There are {data_train.isnull().any().sum()} columns in train dataset with missing values.')

There are 0 columns in train dataset with missing values.

训练集没有缺失值，不需要进行缺失值填充。
同样地查看测试集的基本情况：

print('The shape of test data:', data_test.shape)
data_test.head()

The shape of test data: (782378, 6)

查看测试集各个特征的基本统计信息：

data_test.describe([0.01, 0.025, 0.05, 0.5, 0.75, 0.9, 0.99])

速度特征同样存在异常值。
查看测试集是否存在缺失值：

# 查看测试集缺失值
print(f'There are {data_test.isnull().any().sum()} columns in test dataset with missing values.')

测试集也没有缺失值。

小结

通过对数据集的整体了解，我们发现数据集的特征构成并不复杂，并且都是数值型数据，也没有缺失值，也就是说原始的数据不需要我们进行过多的处理，只需滤掉一些异常值。对于这样一个数据集，我们工作的重点就在于通过分析数据的分布挖掘不同类型渔船的轨迹特点来找到构建新特征的思路，并且通过合适的方式来利用轨迹这一时序信息。
接下来我们就要对数据的分布进行一些可视化的分析。

2. 可视化分析

2.1 渔船轨迹可视化

我们先把三类不同的数据存入不同的文件中：

# 把三类不同数据放到不同的文件中
def get_diff_data():
    Path = '../input/wisdomoceans/total_data.pkl'
    with open(Path, 'rb') as f:
        total_data = pickle.load(f)
    load_save = Load_Save_Data()
    kind_data = ['拖网', '围网', '刺网']
    filenames = ['tuowang_data.pkl', 'weiwang_data.pkl', 'ciwang_data.pkl']
    for i, datax in enumerate(kind_data):
        data_type = [data for data in total_data if data['type'].unique()[0] == datax]
        load_save.save_data(data_type, './' + filenames[i])

get_diff_data()

再分别从三类数据文件中随机取三条渔船可视化他们的轨迹：

# 分别从三类数据文件中，随机读取三条渔船的数据
def get_random_three_traj(type=None):
    np.random.seed(10)
    path = '../input/wisdomoceans/'
    with open (path + type + '.pkl', 'rb') as f:
        data = pickle.load(f)
    data_arrange = np.arange(len(data)).tolist()
    index = random.sample(data_arrange, 3)
    return data[index[0]], data[index[1]], data[index[2]]

# 渔船轨迹可视化
def visualize_three_traj():
    fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(20, 15))
    plt.subplots_adjust(wspace=0.2, hspace=0.2)
    labels = ['tuowang', 'weiwang', 'ciwang']
    for i, filetype in tqdm(enumerate(['tuowang_data', 'weiwang_data', 'ciwang_data'])):
        data1, data2, data3 = get_random_three_traj(filetype)
        for j, datax in enumerate([data1, data2, data3]):
            x_data = datax['x'].loc[-1:].values
            y_data = datax['y'].loc[-1:].values
            axes[i][j-1].scatter(x_data[0], y_data[0], label='start', c='red', s=10, marker='8')
            axes[i][j-1].plot(x_data, y_data, label=labels[i])
            axes[i][j-1].scatter(x_data[len(x_data)-1], y_data[len(y_data)-1], label='end', c='green', s=10, marker='v')
            
            axes[i][j-1].grid(alpha=.2)
            axes[i][j-1].legend(loc='best')
            
    plt.show()

visualize_three_traj()

得到三类数据中随机三条渔船的轨迹：

可以看到：拖网的轨迹通常是从一点到另一点；围网的轨迹通常会形成一个圈；刺网的轨迹通常由多条直线段构成。
但同时也注意到，渔船的轨迹并非是完全按照我们上述的规律走的，例如上图中很多轨迹是较为杂乱的。这是因为渔船并非全程都在执行作业，而是存在停靠点和驻留区域，如下图所示：

PS：上图来自TOP3大佬开源的决赛PPT
因此，如何区分开渔船的机动阶段和POI点、ROI区域，也是特征工程中非常重要的一项工作。

2.2 渔船坐标序列可视化

不同作业类型的渔船的坐标变化是不同的，我们随机取一条渔船来看看它轨迹的x轴和y轴坐标变化情况：

# 随机读取某类数据中一条渔船的数据
def get_random_one_traj(type=None):
    np.random.seed(10)
    path = '../input/wisdomoceans/'
    with open (path + type + '.pkl', 'rb') as f:
        data = pickle.load(f)
    index = np.random.choice(len(data))
    return data[index]

# 可视化x、y变化
def visualize_one_traj_x_y():
    fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(10, 8))
    plt.subplots_adjust(wspace=0.5, hspace=0.5)
    
    # 获取围网数据中某条渔船的数据
    data = get_random_one_traj('weiwang_data')
    x = data['x'].loc[-1:]
    x /= 1000
    
    y = data['y'].loc[-1:]
    y /= 1000
    
    arr1 = np.arange(len(x))
    arr2 = np.arange(len(y))
    
    axes[0].plot(arr1, x, label='x')
    axes[1].plot(arr2, y, label='y')
    axes[0].grid(alpha=0.5)
    axes[0].legend(loc='best')
    axes[1].grid(alpha=0.5)
    axes[1].legend(loc='best')
    plt.show()

visualize_one_traj_x_y()

得到x、y坐标的变化情况：

可以看到，这条围网渔船的x和y轴在同一时刻开始变化，到某一时刻又同时不再改变。并且过程中也存在几段x、y坐标都不变化的时间段，说明渔船停靠在港口中，即存在POI点。

2.3 速度和方向可视化

我们从三类数据文件中各随机取一条渔船的数据，观察它们的速度和方向变化：

def visualize_three_traj_speed_direction():
    fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(20, 15))
    plt.subplots_adjust(wspace=0.3, hspace=0.3)
    
    filetypes = ['tuowang_data', 'weiwang_data', 'ciwang_data']
    speedtypes = ['tuowang_speed', 'weiwang_speed', 'ciwang_speed']
    directiontypes = ['tuowang_direction', 'weiwang_direction', 'ciwang_direction']
    colors = ['pink', 'lightblue', 'lightgreen']
    for i, filename in tqdm(enumerate(filetypes)):
        datax = get_random_one_traj(filename)
        x_data = datax['速度'].loc[-1:].values
        y_data = datax['方向'].loc[-1:].values
        axes[i][0].plot(range(len(x_data)), x_data, label=speedtypes[i], color=colors[i])
        axes[i][0].grid(alpha=0.5)
        axes[i][0].legend(loc='best')
        axes[i][1].plot(range(len(y_data)), y_data, label=directiontypes[i], color=colors[i])
        axes[i][1].grid(alpha=0.5)
        axes[i][1].legend(loc='best')
    plt.show()

visualize_three_traj_speed_direction()

得到三类渔船的速度和方向变化：

可以看到：拖网渔船的速度和方向变化是阶段式的，在一段时间内速度波动大、方向波动小，另一段时间内速度波动小，而方向波动大；围网渔船的速度变化是脉冲式的，间歇性地在某些时间点剧烈变化，而方向在整个过程中都不断变化；刺网渔船的速度和变化是阶梯式的，速度和方向在同一时间段保持不变，然后又同时开始改变。

2.4 速度和方向数据分布可视化

绘制各类渔船的速度和方向数据核密度估计图：

# 绘制训练数据的速度和方向分布图
plt.subplots(nrows=1, ncols=2, figsize=(15, 6))
plt.subplots_adjust(wspace=0.3, hspace=0.5)
types = ['拖网', '围网', '刺网']
labels = ['target==tw', 'target==ww', 'target==cw']
colors = ['red', 'green', 'blue']
for i, t in enumerate(types):
    type_df = data_train[data_train['type']==t]
    plt.subplot(1, 2, 1)
    ax1 = sns.kdeplot(type_df['speed'].values, color=colors[i], shade=True)
    plt.subplot(1, 2, 2)
    ax2 = sns.kdeplot(type_df['direction'].values, color=colors[i], shade=True)
    ax1.legend(labels)
    ax1.set_xlabel('Speed')
    ax2.legend(labels)
    ax2.set_xlabel('Direction')
plt.show()

绘制各类渔船的速度和方向箱线图：

# 使用箱线图进行可视化
def plot_speed_direction2_ditribution():
    fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15, 6))
    plt.subplots_adjust(wspace=0.3, hspace=0.5)
    colors = ['pink', 'lightblue', 'lightgreen']
    
    bplot1 = axes[0].boxplot([df_speeds[0]['speed'].values, df_speeds[1]['speed'].values, df_speeds[2]['speed'].values],
                             vert=True, patch_artist=True, labels=['tw', 'ww', 'cw'])
    bplot2 = axes[1].boxplot([df_directions[0]['direction'].values, df_directions[1]['direction'].values, df_directions[2]['direction'].values],
                             vert=True, patch_artist=True, labels=['tw', 'ww', 'cw'])
    
    for bplot in (bplot1, bplot2):
        for patch, color in zip(bplot['boxes'], colors):
            patch.set_facecolor(color)
    
    axes[0].set_title('Speed')
    axes[1].set_title('Direction')
    plt.show()

plot_speed_direction2_ditribution()

从核密度估计图和箱线图中都可以看出，三类渔船的速度和方向分布存在较明显的区别，因此我们可以利用速度和方向特征来构建新的特征。
同时也注意到，速度和方向的分布都存在离群点，这说明速度和方向中都存在异常值，因此我们需要对异常值进行处理之后，再次分析数据的分布。

3. 异常值处理

在处理异常值之前，为了使可视化的数据更加直观，我们先将x和y轴坐标转换成经纬度：

from pyproj import Proj

# 将xy坐标转换成经纬度
def transform_xy2lonlat(df):
    x = df['lat'].values
    y = df['lon'].values
    p=Proj('+proj=lcc +lat_1=33.88333333333333 +lat_2=32.78333333333333 +lat_0=32.16666666666666 +lon_0=-116.25 +x_0=2000000.0001016 +y_0=500000.0001016001 +datum=NAD83 +units=us-ft +no_defs ')
    df['lon'], df['lat'] = p(y, x, inverse=True)
    return df

new_data_train = data_train.copy()
new_data_train.columns = ['ID', 'lat', 'lon', 'speed', 'direction', 'time', 'type']
new_data_train = transform_xy2lonlat(new_data_train)

new_data_test = data_test.copy()
new_data_test.columns = ['ID', 'lat', 'lon', 'speed', 'direction', 'time']
new_data_test = transform_xy2lonlat(new_data_test)

采用3-sigma算法来判断异常值，但是我们并不把所有异常值都直接剔除掉，对于速度和方向数据，我们可以先将异常值置为空值，再采用多项式插值来进行填充。

# 3-sigma算法判定异常值
def three_sigma(data):
    data_mean = np.mean(data)
    data_std = np.std(data)
    low = data_mean - 3 * data_std
    high = data_mean + 3 * data_std
    judge = []
    for d in data:
        if d < low or d > high:
            judge.append(True)
        else:
            judge.append(False)
    return judge

# 将异常值置为空值
def assign_traj_anomaly_points_nan(df):
    # 速度异常值用空值填充
    is_speed_anomaly = three_sigma(df['speed'])
    df['speed'][is_speed_anomaly] = np.nan
    
    # 方向异常值用空值填充
    is_direction_anomaly = three_sigma(df['direction'])
    df['direction'][is_direction_anomaly] = np.nan
    
    # 纬度和经度异常值直接剔除
    is_lat_anomaly = three_sigma(df['lat'])
    is_lon_anomaly = three_sigma(df['lon'])
    lat_lon_anomaly = np.array([a | b for a, b in zip(is_lat_anomaly, is_lon_anomaly)])
    df = df[~lat_lon_anomaly].reset_index(drop=True)
    
    # 统计异常值个数
    anomaly_cnt = len(is_speed_anomaly) - len(df)
    
    return df, anomaly_cnt

# 判断每个渔船ID轨迹中的异常值，采用多项式插值来填充
train_ids = list(new_data_train['ID'].unique())
train_new = []
train_anomaly_cnts = []
for ID in tqdm(train_ids):
    id_df = new_data_train[new_data_train['ID']==ID]
    id_df, anomaly_cnt = assign_traj_anomaly_points_nan(id_df)
    # 对速度和方向异常值进行二阶插值
    id_df = id_df.interpolate(method='polynomial', axis=0, order=2)
    # 填充剩下的空值
    id_df = id_df.fillna(method='bfill')
    id_df = id_df.fillna(method='ffill')
    # 设置阈值过滤掉不合理的值
    id_df['speed'] = id_df['speed'].clip(0, 23)
    id_df['direction'] = id_df['direction'].clip(0, 360)
    # 统计每个id的异常值个数
    train_anomaly_cnts.append(anomaly_cnt)
    train_new.append(id_df)
new_data_train = pd.concat(train_new)

我们选取ID为3265的渔船来看一下异常值过滤前后的轨迹，这是一条拖网作业的渔船：

# 单条渔船轨迹可视化
def visualize_one_traj(df, new_df, ID):
    fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(20, 8))
    plt.subplots_adjust(wspace=0.2, hspace=0.2)
    id_data, new_id_data = df[df['ID']==ID], new_df[new_df['ID']==ID]
    labels = ['before', 'later']
    for i, datax in enumerate([id_data, new_id_data]):
        x_data = datax['lat'].loc[-1:].values
        y_data = datax['lon'].loc[-1:].values
        axes[i].scatter(x_data[0], y_data[0], label='start', c='red', s=10, marker='8')
        axes[i].plot(x_data, y_data, label=labels[i])
        axes[i].scatter(x_data[len(x_data)-1], y_data[len(y_data)-1], label='end', c='green', s=10, marker='v')
            
        axes[i].grid(alpha=0.2)
        font = {'family': 'Times New Roman', 'weight': 'normal', 'size': 20}
        axes[i].legend(loc='best', prop=font)
            
    plt.show()

visualize_one_traj(data_train, new_data_train, 3265)

轨迹中最底下位置跨度大的一部分被过滤掉了，过滤后的轨迹是点到点的，可以很明显地看出是拖网作业。

4. 滤噪后的数据可视化分析

处理了异常值后，我们重新对速度和方向分布进行可视化分析，得到的核密度估计图和箱线图如下：


滤噪后的速度和方向分布没有了离群点，可以明显地看到三类渔船数据分布是不同的：

• 刺网渔船速度低，方向角小
• 拖网渔船速度高，方向角小
• 围网渔船速度均值比刺网渔船高，方向角大

同样绘制不同类型以及总数据的经纬度分布核密度估计图：

上图表明，不同类型渔船的位置分布也存在明显的差异性，可以由此构建相关的特征。

5. 总结

经过以上的数据分析，我们可以得出以下结论：

• 数据中存在异常值，可通过3-sigma算法和多项式插值进行处理
• 渔船轨迹存在POI点和ROI区域，需要对这些信息进行挖掘
• 可根据速度、方向以及经纬度这些有区分度的特征构造新的特征