AIS数据下载并处理（python）

最近由于项目需求，需要下载AIS数据，其实也就是船舶位置数据。

首先通过百度发现目前一些网站下载的AIS数据都是收费的，如船达通等等。

船达通是只能查找近一个月的船舶轨迹信息，下载需要支付收费。

但是最后我终于找到了一个免费下载AIS数据的网站，并且下面附上了处理说明，可以通过python代码批量处理

通过以下的文档可以详细了解如何下载数据，将AIS数据中的时间、经度、纬度、航速提取出来。

最后通过ArcMap处理成为shp，方便可视化。

我是从中国海洋卫星数据网站下载的。

网站地址：https://osdds.nsoas.org.cn/#/

需要注册后使用。

有关网站注册和注册登录不上的问题，可以去邮件网络官方，我当时用没问题，不了解。

另外不方便借账号。

一、数据说明文档

首先从首页的下载AIS 0-1级产品数据格式说明V1.0.pdf

打开后可以看到，HY-C/D的L1A产品是可以直接使用的，因此从网站上下载该产品。

二、数据下载

点击网站首页的数据获取，点击左边的海洋水色卫星数据获取可以得到右边这个界面。

按需求选择。

加入订单，下载即可，我觉得有点麻烦的是，批量勾选，每勾选一次自动跳回第一页= =。

最后就等等审批（应该很快，几分钟？），然后就可以下载了/

推荐使用DownThemAll批量下载器，火狐的插件。

三、数据处理

数据下载下来，压缩包格式，解压后，L1A后缀文件才是AIS数据，是csv格式的，并且用逗号分隔，可以用excel打开。

那么批量处理我主要是分成两个程序了。

1、解压

# -*- coding: utf-8 -*-
"""
Created on Sat Feb  6 22:02:47 2021
解压tar.gz文件
@author: 就是一只白
"""

import tarfile
import os
import glob
def untar(fname, dirs):
    t = tarfile.open(fname)
    t.extractall(path = dirs) 
rootdir = "H:\\AIS\\20210104_20210119\\"
outfolder='H:\\AIS\\20210104_20210119\\process\\'
files = glob.glob(rootdir+"*.tar.gz")
#filelist = os.listdir(rootdir)
for file in files:
    print(file)
    untar(file,outfolder)

2、提取时间、位置等数据

报文的解析可以看pdf文件，CSV文件每行代表一个报文消息，每一列是该条信息的一个字段（如速度、时间）。一共有27个报文类型，不同类型的报文后面的信息字段不同。

我下面的代码只筛选了几种进行处理，因为通过excel表筛选发现其他的我这边没有。

以下代码可以获取两个csv文件，代码注释写了两个文件的区别。

最后其实只需要用position的文件即可，数据有四项：时间、航速、经度、纬度

通过修改下方的mmsi号修改查找的船舶

# -*- coding: utf-8 -*-
"""
Created on Sun Feb  7 10:24:29 2021

@author: jiushiyizhibai
"""


import os
import glob
from pandas.core.frame import DataFrame
import csv
import numpy as np
import datetime
#按照mmsi号查找报文信息
rootdir = "H:\\AIS\\20210104_20210119\\process\\"
outfolder='H:\\AIS\\20210104_20210119\\csv\\'
mmsi='309351666'
files = glob.glob(rootdir+"*.l1a")
result=[]
#filelist = os.listdir(rootdir)
for file in files:
    print(file)
    table=open(file,'r')
    s = [i[:-1].split(',') for i in table.readlines()]
    s=[n for n in s if n[5]==mmsi]
    print(len(s))
    if len(s)>0 :
        for i in s:
            result.append(i)
data=DataFrame(result)
time=[]
speed=[]
lon=[]
lat=[]
for row in data.itertuples():
    if row[4]=='1' or row[4]=='3':
        time.append(datetime.datetime.strptime(row[1],'%Y-%m-%d %H:%M:%S'))
        speed.append(float(row[9]))
        lon.append(float(row[11]))
        lat.append(float(row[12]))
    elif row[4]=='27':
        time.append(datetime.datetime.strptime(row[1],'%Y-%m-%d %H:%M:%S'))
        speed.append(float(row[12]))
        lon.append(float(row[10]))
        lat.append(float(row[11]))
    else:
        print("没有查到此类报文：",row[4])
        
presult = DataFrame({'time':time,
                      'speed':speed,
                      'lon':lon,
                      'lat':lat})
presult=presult.sort_values(by='time')
#['time','speed','lon','lat']
#print(data)
data=data.values
#输出文件csv，这里是只输出时间速度和位置
presult.to_csv(outfolder+mmsi+'_AIS_position.csv',index = False)
#输出所有报文信息的文件
csvfile = open(outfolder+mmsi+'_AIS.csv', 'w')  #打开方式还可以使用file对象
writer = csv.writer(csvfile)
writer.writerows(data)
csvfile.close()

3、将csv转成点和线的shapefile轨迹文件（ArcMap)

如何将Excel经纬度坐标数据导入ArcGIS

看上面的即可，点转线用的是点集转线这个功能。

PS:XY建立点数据的时候，选择坐标系统，不要选择投影系统！！不然就会出现点不知在哪的情况