1.在NOAA(national oceanic and atmospheric administration)官网下载迈阿密港口船舶进出数据,时间为2008/12/31 23:59:00到2009/1/4 23:59:00,时间间隔为1分钟。以下都是轨迹点数据:
2.打开属性表,最主要的字段为时间BaseDateTime,船舶号MMSI(相当于车辆的id编号)
3.根据船舶号MMSI,利用轨迹追踪工具脚本,根据时间顺序生成轨迹线数据
代码也是由NOAA提供的,船舶轨迹点数据生成轨迹线数据
"""
August 21, 2015
Marine Cadastre
www.marinecadastre.gov
NOAA Office for Coastal Management
1234 South Hobson Avenue
Charleston, SC 29405
843-740-1200
ocm.mmc@noaa.gov
Software design and development by
RPS ASA, Inc.
55 Village Square Drive
Wakefield, RI 02879
IMSG
3206 Tower Oaks Blvd
Rockville, MD 20852
"""
########################################################################################################
#########################################################################################################
#########################################################################################################
## Input Paramaters
#Path to input points (broadcast feature class)
sInputFile = r"D:\Miami.gdb\Broadcast"
#ID field name (can keep as is most likely)
sID_Field = "MMSI"
#DateTime field name (can keep as is most likely)
sDT_Field = "BaseDateTime"
#Path to output workspace (file geodatabase)
sOutWorkspace = r"D:\Miami.gdb"
#name of output trackline feature class
sOutName = "TrackLine"
#Break Track line method option (comment/uncomment the selected option)
sLineMethod = "Maximum Time and Distance"
#sLineMethod = "Time Interval"
#Maximum Time in minutes (keep as string in quotes, use 0 if using Time Interval method, default = 60)
sMaxTimeMinutes = "60"
#Maximum Distance in miles (keep as string in quotes, use 0 if using Time Interval method, default = 6)
sMaxDistMiles = "6"
#Time Interval in hours (keep as string in quotes, use 0 if using Maximum Time and Distance method, default = 24)
sIntervalHours = "0"
#Include Voyage Attributes (comment/uncomment the selected option)
#sAddVoyageData = "true"
sAddVoyageData = "false"
#Voyage Table (path to voyage table) - OPTIONAL
voyageTable = ""
#Voyage Matching Option (comment/uncomment the selected option) - OPTIONAL
#voyageMethod = "Most Frequent"
voyageMethod = "First"
#Include Vessel Attributes (comment/uncomment the selected option)
#sAddVesselData = "true"
sAddVesselData = "false"
#Vessel Table (path to vessel table) - OPTIONAL
vesselTable = ""
#########################################################################################################
#########################################################################################################
#########################################################################################################
import arcpy, os, sys, traceback, datetime, math, re
import multiprocessing
import time as t
from functools import partial
iMaxTimeMinutes = 0
iMaxDistMiles = 0
iIntervalHours = 0
spatialRef_CURSOR = None
#########################################################################################################
class Timer:
def __init__(self):
self.startTime = t.time()
def End(self):
self.endTime = t.time()
self.elapTime = self.endTime - self.startTime
return self.elapTime
class Progress:
def __init__(self, total):
self.iTotal = float(total)
self.prog = 0
sys.stdout.write('\r 0%')
def Update(self, current):
perc = round(current/self.iTotal * 100.0,0)
if perc > self.prog:
self.prog = perc
sys.stdout.write('\r '+str(perc)+"%")
sys.stdout.flush()
def CheckAISDatabase():
"""Performs a series of checks to see if the source points are a true AIS database with vessel and voyage data."""
bAISdb = True
desc = arcpy.Describe(sInputFile)
input_wkspace = desc.Path
sInName = desc.Name
if os.path.splitext(input_wkspace)[1] <> ".gdb":
bAISdb = False
else:
if sInName.find("Broadcast") < 0:
bAISdb = False
else:
sVesselName = re.sub("Broadcast", "Vessel", sInName)
sVoyageName = re.sub("Broadcast", "Voyage", sInName)
if arcpy.Exists(input_wkspace+"\\"+sVesselName) == False:
bAISdb = False
else:
if arcpy.Exists(input_wkspace+"\\"+sVoyageName) == False:
bAISdb = False
else:
if arcpy.Exists(input_wkspace+"\\BroadcastHasVessel") == False:
bAISdb = False
else:
if arcpy.Exists(input_wkspace+"\\BroadcastHasVoyage") == False:
bAISdb = False
return bAISdb
def FormatCounts(iInCount):
"""Formats various counts to strings with comma separators, for use in dialog messages."""
sOutCount = re.sub("(\d)(?=(\d{3})+(?!\d))", r"\1,", "%d" % iInCount)
return sOutCount
def GetVersion():
"""Gets the version of ArcGIS being used. Faster methods to read input datasets are used when 10.1 is the version in use."""
dicInstallInfo = arcpy.GetInstallInfo()
return dicInstallInfo["Version"]
def CreateDataDict_NEW():
"""Process to read through the input point feature class. Point information is stored in a Python dictionary in memory.
All point attributes and coordinates (expect date/time) are read into Numpy Array.
The date information is read using the Data Access module search cursor.
The two parts are then combined based on the OID into the Python dictionary"""
print "\n Reading input point data..."
iTotPoints = int(arcpy.GetCount_management(sInputFile).getOutput(0))
timer = Timer()
prog = Progress(iTotPoints*3) #once for array, date dict, merge
tally = 0
i=0
if bRealAISdb == True:
fields = ["SHAPE@XY",sID_Field,"VoyageID","OID@"]
fields2 = ["OID@",sDT_Field]
##FILTERING MMSI CODES THAT ARE ZERO
where = sID_Field+" > 0"
else:
fields = ["SHAPE@XY",sID_Field,"OID@"]
fields2 = ["OID@",sDT_Field]
where = None
#Store all data (except date) into numpy array
ptData_Array = arcpy.da.FeatureClassToNumPyArray(sInputFile,fields,where,spatialRef_CURSOR)
tally+=iTotPoints
prog.Update(tally)
#Get date using search cursor, store in temporary dictionary
dateDict = {}
with arcpy.da.SearchCursor(sInputFile,fields2,where,spatialRef_CURSOR,False) as rows:
for row in rows:
sOID = row[0]
dDateTime = row[1]
dateDict[sOID] = dDateTime
tally+=1
prog.Update(tally)
#to account for points without MMSI
tally=iTotPoints*2
prog.Update(tally)
#combine array and date dictionary into one final dictionary
dataDict = {}
for i in xrange(0,len(ptData_Array[sID_Field])):
if bRealAISdb == True:
sID, oid, voyageID, geo = ptData_Array[sID_Field][i], ptData_Array["OID@"][i], ptData_Array["VoyageID"][i], ptData_Array["SHAPE@XY"][i]
else:
sID, oid, geo, voyageID = ptData_Array[sID_Field][i], ptData_Array["OID@"][i], ptData_Array["SHAPE@XY"][i], None
if dataDict.has_key(sID):
dataDict[sID][dateDict[oid]] = [geo[0],geo[1],voyageID]
else:
dataDict[sID] = {dateDict[oid]:[geo[0],geo[1],voyageID]}
tally+=1
prog.Update(tally)
del sID, oid, geo, voyageID
del dateDict, ptData_Array
print "\r Read "+FormatCounts(iTotPoints)+" points"
print " Read Complete:" + str(timer.End())
del timer, prog
return dataDict, iTotPoints
def SplitDataDictEqually(input_dict, parts, total_points):
"""Splits the point data dictionary into equal parts, for separate processing.&