Datawhale 智慧海洋建设-Task2 数据分析

此部分为智慧海洋建设竞赛的数据分析模块，通过数据分析，可以熟悉数据，为后面的特征工程做准备，欢迎大家后续多多交流。

赛题：智慧海洋建设 - `$ V0 S7 P: _% [* T, v

数据分析的目的:

EDA的主要价值在于熟悉整个数据集的基本情况(缺失值、异常值)，来确定所获得数据集可以用于接下来的机器学习或者深度学习使用。了解特征之间的相关性、分布，以及特征与预测值之间的关系。为进行特征工程提供理论依据。

项目地址：https://github.com/datawhalechina/team-learning-data-mining/tree/master/wisdomOcean比赛地址：https://tianchi.aliyun.com/competition/entrance/231768/introduction?spm=5176.12281957.1004.8.4ac63eafE1rwsY

2.1 学习目标

学习如何对数据集整体概况进行分析，包括数据集的基本情况(缺失值、异常值)学习了解变量之间的相互关系、变量与预测值之间的存在关系。完成相应学习打卡任务

2.2 内容介绍

数据总体了解读取数据集并了解数据集的大小，原始特征维度；通过info了解数据类型；粗略查看数据集中各特征的基本统计量缺失值和唯一值查看数据缺失值情况查看唯一值情况

数据特性和特征分布

# h3 q! F0 K+ D8 b0 K7 K

三类渔船轨迹的可视化坐标序列可视化三类渔船速度和方向序列可视化三类渔船速度和方向的数据分布

作业一：剔除异常点后画图

import pandas as pd

import geopandas as gpd

9 V4 P" x$ ]/ p: j+ ] s2 l

from pyproj import Proj

7 k# V; ^+ c; R" ~* I" H

from keplergl import KeplerGl

8 Z1 D9 A' o S/ w5 z' x J+ R4 F

from tqdm import tqdm

+ T5 k2 v) z6 U9 r

import os

! x* P* @6 w7 Q# |3 r0 p8 R0 o3 j

import matplotlib.pyplot as plt

( ]5 B1 H; d" B

import shapely

6 m- o3 C$ N2 {' f; [' H: t

import numpy as np

! U) r- Q( d/ n/ B4 k

from datetime import datetime

9 u, K" e; I* Q8 L* t: I3 u

import warnings

& n4 f1 O7 r5 C: j

warnings.filterwarnings(ignore)

' r o& B% a0 V2 k

plt.rcParams[font.sans-serif] = [SimSun] # 指定默认字体为新宋体。

! \# Y3 h0 u6 d3 T2 }( m5 ~

plt.rcParams[axes.unicode_minus] = False # 解决保存图像时 负号- 显示为方块和报错的问题。

#获取文件夹中的数据

def get_data(file_path,model):

# M2 O* \$ e5 U

assert model in [train, test], {} Not Support this type of file.format(model)

# J6 T3 H) L& _0 ` e5 h2 z

paths = os.listdir(file_path)

# print(len(paths))

+ y- i7 P$ E# b2 s

tmp = []

2 z! o7 X6 c0 k1 u3 |, O: W

for t in tqdm(range(len(paths))):

! g2 E& u: A$ q

p = paths[t]

with open({}/{}.format(file_path, p), encoding=utf-8) as f:

+ Q% I Z; ]$ k7 [' @6 ?

next(f)

for line in f.readlines():

tmp.append(line.strip().split(,))

tmp_df = pd.DataFrame(tmp)

if model == train:

tmp_df.columns = [ID, lat, lon, speed, direction, time, type]

+ l% H& \) b, Q

else:

4 D; B1 q# p. t1 l

tmp_df[type] = unknown

tmp_df.columns = [ID, lat, lon, speed, direction, time, type]

% ^: V6 z- g+ c( m; U. ^

tmp_df[lat] = tmp_df[lat].astype(float)

tmp_df[lon] = tmp_df[lon].astype(float)

tmp_df[speed] = tmp_df[speed].astype(float)

! ?0 M# S: y7 t0 S& @' R2 [$ y

tmp_df[direction] = tmp_df[direction].astype(int)#如果该行代码运行失败，请尝试更新pandas的版本

return tmp_df

5 P. F( f" D$ A5 O! D& M3 b7 _* U, o

# 平面坐标转经纬度，供初赛数据使用

) V- W \; U% X0 `( c5 M# M/ w

# 选择标准为NAD83 / California zone 6 (ftUS) (EPSG:2230)，查询链接：CS2CS - Transform Coordinates On-line - MyGeodata Cloud

def transform_xy2lonlat(df):

8 I; i Y' b9 d; v6 e0 z

x = df[lat].values

$ _; L& V( D! W, Q

y = df[lon].values

; C- t5 p( d- j8 x* h

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 )

" w: j7 {/ r4 ], a' R* _" |

df[lon], df[lat] = p(y, x, inverse=True)

2 R4 x+ a+ o: N4 h

return df

/ k1 [) [3 n1 ]6 ~ q2 D" M1 n# A

#修改数据的时间格式

2 V" ]1 I: c' u0 _* L% l

def reformat_strtime(time_str=None, START_YEAR="2019"):

( f0 b4 r0 s/ C

"""Reformat the strtime with the form 08 14 to START_YEAR-08-14 """

1 z1 F! u( f' n9 c+ e+ a

time_str_split = time_str.split(" ")

time_str_reformat = START_YEAR + "-" + time_str_split[0][:2] + "-" + time_str_split[0][2:4]

# Z) y4 A* j" _

time_str_reformat = time_str_reformat + " " + time_str_split[1]

: P! \, d2 g5 q& \0 l. v# N6 T5 ^3 a+ B

# time_reformat=datetime.strptime(time_str_reformat,%Y-%m-%d %H:%M:%S)

return time_str_reformat

6 ~3 n" Y4 f" w+ L4 x, z# G+ s0 T

#计算两个点的距离

def haversine_np(lon1, lat1, lon2, lat2):

8 T+ e- \1 G, \' {: j# @' M

lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

+ N! _+ m- F3 c

dlon = lon2 - lon1

& s7 ~ P9 l$ w/ }

dlat = lat2 - lat1

/ L4 B3 q; Z# Z) \8 f8 B M

a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2

c = 2 * np.arcsin(np.sqrt(a))

7 k& y7 D5 G4 p0 f, \; y1 k

km = 6367 * c

return km * 1000

def compute_traj_diff_time_distance(traj=None):

"""Compute the sampling time and the coordinate distance."""

# 计算时间的差值

) R( g) L: G3 X& Y3 M$ A" A

time_diff_array = (traj["time"].iloc[1:].reset_index(drop=True) - traj[

; a7 B1 I8 q% p

"time"].iloc[:-1].reset_index(drop=True)).dt.total_seconds() / 60

4 A4 t4 A. g9 i/ `: S7 e! S

# 计算坐标之间的距离

dist_diff_array = haversine_np(traj["lon"].values[1:], # lon_0

& B0 o; N3 P1 b$ p: R

traj["lat"].values[1:], # lat_0

traj["lon"].values[:-1], # lon_1

traj["lat"].values[:-1] # lat_1

)

# 填充第一个值

time_diff_array = [time_diff_array.mean()] + time_diff_array.tolist()

& r3 v- E7 ?, w0 i; U

dist_diff_array = [dist_diff_array.mean()] + dist_diff_array.tolist()

traj.loc[list(traj.index),time_array] = time_diff_array

traj.loc[list(traj.index),dist_array] = dist_diff_array

1 g" `. m+ Q9 L [; g& ^+ m

return traj

#对轨迹进行异常点的剔除

: w% N' H$ o0 T3 f7 o& q9 U/ P7 S6 J

def assign_traj_anomaly_points_nan(traj=None, speed_maximum=23,

( c# {4 r$ @6 F C

time_interval_maximum=200,

coord_speed_maximum=700):

"""Assign the anomaly points in traj to np.nan."""

def thigma_data(data_y,n):

' B1 a2 X2 d) E

data_x =[i for i in range(len(data_y))]

ymean = np.mean(data_y)

1 x0 T4 G2 y, j4 D; G

ystd = np.std(data_y)

) x2 C; A$ D! G" ~5 i; n1 w

threshold1 = ymean - n * ystd

threshold2 = ymean + n * ystd

judge=[]

" ^- B! x1 G! X% Y& o! N7 V$ \

for data in data_y:

if (data < threshold1)|(data> threshold2):

judge.append(True)

7 Q. s# b7 |6 V' I/ O% H! | T# d. m

else:

( Y; o2 J% W$ ~5 }% J6 L& @3 j

judge.append(False)

return judge

# Step 1: The speed anomaly repairing

is_speed_anomaly = (traj["speed"] > speed_maximum) | (traj["speed"] < 0)

7 g% B6 m8 R8 G$ G7 |$ q

traj["speed"][is_speed_anomaly] = np.nan

) S$ @9 U" H# m0 {$ ?- c/ X

# Step 2: 根据距离和时间计算速度

( L$ x1 o# P) z2 ]; q. o

is_anomaly = np.array([False] * len(traj))

* D5 e8 I4 Q$ Q8 |. T" \$ l

traj["coord_speed"] = traj["dist_array"] / traj["time_array"]

9 M$ k4 }5 q4 c& J( }# N

# Condition 1: 根据3-sigma算法剔除coord speed以及较大时间间隔的点

is_anomaly_tmp = pd.Series(thigma_data(traj["time_array"],3)) | pd.Series(thigma_data(traj["coord_speed"],3))

is_anomaly = is_anomaly | is_anomaly_tmp

is_anomaly.index=traj.index

1 e$ U E1 E( \% @, |

# Condition 2: 轨迹点的3-sigma异常处理

traj = traj[~is_anomaly].reset_index(drop=True)

# d9 C% F* N- A

is_anomaly = np.array([False] * len(traj))

2 q2 i8 g# y- Q; b/ H

if len(traj) != 0:

( f+ ^- R" H; U

lon_std, lon_mean = traj["lon"].std(), traj["lon"].mean()

lat_std, lat_mean = traj["lat"].std(), traj["lat"].mean()

lon_low, lon_high = lon_mean - 3 * lon_std, lon_mean + 3 * lon_std

lat_low, lat_high = lat_mean - 3 * lat_std, lat_mean + 3 * lat_std

is_anomaly = is_anomaly | (traj["lon"] > lon_high) | ((traj["lon"] < lon_low))

is_anomaly = is_anomaly | (traj["lat"] > lat_high) | ((traj["lat"] < lat_low))

traj = traj[~is_anomaly].reset_index(drop=True)

/ `/ T x: Y$ I: N( `

return traj, [len(is_speed_anomaly) - len(traj)]

* D; K: I; V$ U5 M7 G% A6 h0 H

df=get_data(rC:\Users\admin\hy_round1_train_20200102,train)

3 W9 L. b% e! g$ g3 `/ r+ V0 H( m

#对轨迹进行异常点剔除，对nan值进行线性插值

5 f9 b1 x- h7 Z+ e( r( p" N9 c. k

ID_list=list(pd.DataFrame(df[ID].value_counts()).index)

8 s' X6 Q' v9 r0 F

DF_NEW=[]

?- z1 f( L3 g

Anomaly_count=[]

$ {; w& B* C3 Y% \

for ID in tqdm(ID_list):

df_id=compute_traj_diff_time_distance(df[df[ID]==ID])

& _+ K$ l9 S* o: B+ Q. h1 H# W

df_new,count=assign_traj_anomaly_points_nan(df_id)

- q% T, q6 p5 a" [% M3 c+ Z6 ]

df_new["speed"] = df_new["speed"].interpolate(method="linear", axis=0)

df_new = df_new.fillna(method="bfill")

% C# q O. d8 O0 _- s& ]( P

df_new = df_new.fillna(method="ffill")

! d/ j1 D7 A# c7 q' |* g" i

df_new["speed"] = df_new["speed"].clip(0, 23)

Anomaly_count.append(count)#统计每个id异常点的数量有多少

DF_NEW.append(df_new)

#将数据写入到pkl格式

load_save = Load_Save_Data()

load_save.save_data(DF_NEW,"C:/Users/admin/wisdomOcean/data_tmp1/total_data.pkl")

' V1 P7 e4 [; W( f5 M5 i# `

#### 三类渔船速度和方向可视化

$ c: K7 D- {5 z

# 把训练集的所有数据,根据类别存放到不同的数据文件中

def get_diff_data():

Path = "C:/Users/admin/wisdomOcean/data_tmp1/total_data.pkl"

" E' V$ L4 z! w7 O) V8 m: _: Q+ {

with open(Path,"rb") as f:

total_data = pickle.load(f)

load_save = Load_Save_Data()

4 i! N3 D) Y+ ]. ^; C7 Z; u: a

kind_data = ["刺网","围网","拖网"]

! t0 z+ U4 Z7 I

file_names = ["ciwang_data.pkl","weiwang_data.pkl","tuowang_data.pkl"]

. Q; l7 V; _5 Y

for i,datax in enumerate(kind_data):

data_type = [data for data in total_data if data["type"].unique()[0] == datax]

* l) y+ `, f1 @

load_save.save_data(data_type,"C:/Users/admin/wisdomOcean/data_tmp1/" + file_names[i])

; C( [ v$ c0 r$ ^" N# A3 |

get_diff_data()

#对轨迹进行异常点剔除，对nan值进行线性插值

ID_list=list(pd.DataFrame(df[ID].value_counts()).index)

/ W& j1 w9 Y( \6 R4 v2 p

DF_NEW=[]

$ f0 s" x' y `) X# P8 l9 ]( x! e# p

Anomaly_count=[]

- e! U Q5 t2 l% a

for ID in tqdm(ID_list):

+ H& ]6 {2 U* ], A" T( ]& E0 i" J+ T

df_id=compute_traj_diff_time_distance(df[df[ID]==ID])

; Y1 \) }3 t" o% t8 V- H! \7 d

df_new,count=assign_traj_anomaly_points_nan(df_id)

df_new["speed"] = df_new["speed"].interpolate(method="linear", axis=0)

3 X% F# f8 j4 p! }0 g

df_new = df_new.fillna(method="bfill")

- C4 I3 L& c; B( J) @" [9 `

df_new = df_new.fillna(method="ffill")

df_new["speed"] = df_new["speed"].clip(0, 23)

5 _: P1 J( u2 D: N: {: C

Anomaly_count.append(count)#统计每个id异常点的数量有多少

' g+ w! l: `0 E! \1 G" n$ J/ @0 n

DF_NEW.append(df_new)

* q# ?( G6 ?9 S

# 每类轨迹，随机选取某个渔船，可视化速度序列和方向序列

def visualize_three_traj_speed_direction():

5 C1 O8 ^$ K6 Z! c+ c0 U' t" M( A

fig,axes = plt.subplots(nrows=3,ncols=2,figsize=(20,15))

) i4 k" S9 {: t5 \

plt.subplots_adjust(wspace=0.3,hspace=0.3)

# 随机选出刺网的三条轨迹进行可视化

0 A9 K; z! D; Q4 u) X

file_types = ["ciwang_data","weiwang_data","tuowang_data"]

speed_types = ["ciwang_speed","weiwang_speed","tuowang_speed"]

! @* C# g+ h, W" O! B4 p7 s3 q

doirections = ["ciwang_direction","weiwang_direction","tuowang_direction"]

colors = [pink, lightblue, lightgreen]

6 I: H- V* v" Y$ v

for i,file_name in tqdm(enumerate(file_types)):

datax = get_random_one_traj(type=file_name)

. Z7 n- [: H! _6 @

x_data = datax["速度"].loc[-1:].values

2 P& G: e9 B1 m$ T* g

y_data = datax["方向"].loc[-1:].values

# {" `( J! T# ]( @. ]- G

axes[i][0].plot(range(len(x_data)), x_data, label=speed_types[i], color=colors[i])

/ U) ^, _; x0 S# W. J i) b

axes[i][0].grid(alpha=2)

|3 Z0 T! o+ _: G

axes[i][0].legend(loc="best")

$ Z! l: y* n0 t9 w; F# g

axes[i][1].plot(range(len(y_data)), y_data, label=doirections[i], color=colors[i])

axes[i][1].grid(alpha=2)

$ z- r4 m; C* b

axes[i][1].legend(loc="best")

plt.show()

& w2 o1 Q8 _0 j' l

visualize_three_traj_speed_direction()

. T6 ]# t& c& y

, r% m5 A5 h) F" x) ]% c

作业二：相关性分析。

3 X4 @8 ~! c* z' i

data_train.loc[data_train[type]==刺网,type_id]=1

* g& v' [8 t& ^6 j

data_train.loc[data_train[type]==围网,type_id]=2

* @/ Y& C$ B5 X6 H8 }6 \- a$ C: ~) }

data_train.loc[data_train[type]==拖网,type_id]=3

f, ax = plt.subplots(figsize=(9, 6))

7 M& T# M9 H' }1 h3 A; X

ax = sns.heatmap(np.abs(df.corr()),annot=True)

plt.show()

7 w( b" ^9 w) ?! ?7 p' k

从图中可以清楚看到，经纬度和速度跟类型相关性比较大。

7 a, q% J- ?6 Q9 [. F 0 g. W9 B% o* p5 f : ~/ p8 B2 p u2 i8 L3 k