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."""

return_list = [dict() for idx in xrange(parts)]

min_pts_per_part = math.ceil(total_points/parts)

idx = 0

pts_added = 0

for k,v in input_dict.iteritems():

if pts_added < min_pts_per_part:

return_list[idx][k] = v

pts_added+=len(v.keys())

else:

idx += 1

pts_added = 0

return_list[idx][k] = v

pts_added+=len(v.keys())

return return_list

def CreateOuputDataset(parts=0):

"""Creates an empty feature class to store the track polylines. Adds the required fields, based on the options selected by the user."""

print "\n Building output track line feature class..."

tmp = os.path.join('in_memory', 'template')

#create a feature class in memory, to use as a template for all feature classes created

arcpy.CreateFeatureclass_management('in_memory','template',"POLYLINE","","DISABLED","DISABLED",sInputFile)

arcpy.AddField_management(tmp,sID_Field,"LONG",)

arcpy.AddField_management(tmp,"TrackStartTime","DATE",)

arcpy.AddField_management(tmp,"TrackEndTime","DATE",)

if bRealAISdb == True:

if sAddVoyageData == "true":

arcpy.AddField_management(tmp,"VoyageID","LONG",)

arcpy.AddField_management(tmp,"Destination","TEXT")

arcpy.AddField_management(tmp,"Cargo","LONG",)

arcpy.AddField_management(tmp,"Draught","LONG",)

arcpy.AddField_management(tmp,"ETA","DATE",)

arcpy.AddField_management(tmp,"StartTime","DATE",)

arcpy.AddField_management(tmp,"EndTime","DATE",)

if sAddVesselData == "true":

arcpy.AddField_management(tmp,"IMO","LONG",)

arcpy.AddField_management(tmp,"CallSign","TEXT",)

arcpy.AddField_management(tmp,"Name","TEXT",)

arcpy.AddField_management(tmp,"VesselType","LONG",)

arcpy.AddField_management(tmp,"Length","LONG",)

arcpy.AddField_management(tmp,"Width","LONG",)

arcpy.AddField_management(tmp,"DimensionComponents","TEXT",)

#create a temporary FGDB and feature class for each separate process to avoid locking issues

if parts > 1:

for i in range(1,parts+1):

tmp_wkspc_path = os.path.split(sOutWorkspace)[0]

arcpy.CreateFileGDB_management(tmp_wkspc_path, "temp" + str(i) + ".gdb")

tmp_fgdb = os.path.join(tmp_wkspc_path, "temp" + str(i) + ".gdb")

arcpy.CreateFeatureclass_management(tmp_fgdb,sOutName+str(i),"POLYLINE",tmp,"DISABLED","DISABLED",sInputFile)

#create one feature class since no multiprocessing being used

else:

arcpy.CreateFeatureclass_management(sOutWorkspace,sOutName,"POLYLINE",tmp,"DISABLED","DISABLED",sInputFile)

#delete the temporary template dataset in memory

arcpy.Delete_management(tmp)

del tmp

def GetDistance(prevX,prevY,currX,currY):

"""Calculates the distance between two latitude/longitude coordinate pairs, using the Vincenty formula for ellipsoid based distances.

Returns the distance in miles."""

#Vincenty Formula

#Copyright 2002-2012 Chris Veness

#http://www.movable-type.co.uk/scripts/latlong-vincenty.html

a = 6378137

b = 6356752.314245

f = 1/298.257223563

L = math.radians(currX - prevX)

U1 = math.atan((1-f)*math.tan(math.radians(prevY)))

U2 = math.atan((1-f)*math.tan(math.radians(currY)))

sinU1 = math.sin(U1)

sinU2 = math.sin(U2)

cosU1 = math.cos(U1)

cosU2 = math.cos(U2)

lam = L

lamP = 9999999999

iter_count = 0

while abs(lam-lamP) > 1e-12 and iter_count <= 100:

sinLam = math.sin(lam)

cosLam = math.cos(lam)

sinSigma = math.sqrt((cosU2*sinLam)*(cosU2*sinLam) + (cosU1*sinU2-sinU1*cosU2*cosLam)*(cosU1*sinU2-sinU1*cosU2*cosLam))

if sinSigma == 0:

return 0

cosSigma = sinU1*sinU2+cosU1*cosU2*cosLam

sigma = math.atan2(sinSigma,cosSigma)

sinAlpha = cosU1*cosU2*sinLam/sinSigma

cosSqAlpha = 1 - sinAlpha*sinAlpha

if cosSqAlpha == 0: #catches zero division error if points on equator

cos2SigmaM = 0

else:

cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha

if cos2SigmaM == None:

cos2SigmaM = 0

C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha))

lamP = lam

lam = L+(1-C)*f*sinAlpha*(sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)))

iter_count+=1

uSq = cosSqAlpha*(a*a - b*b)/(b*b)

A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))

B = uSq/1024*(256+uSq*(-128+uSq*(74-47*uSq)))

deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)))

s = b*A*(sigma-deltaSigma)

#convert s in meters to miles

s_miles = s*0.0006213712

return s_miles

def GetTimeDifference(prevDT,currDT):

"""Calculates the difference between to date and time variables. The difference is returned in minutes."""

timeDelta = currDT - prevDT

totMinutes = (timeDelta.days*1440) + (timeDelta.seconds/60)

return totMinutes

def MultiProcess_Main(dataDict, iTotPoints):

"""Main process to split the data and run seperate processes to build the track lines for each available CPU."""

print "\n Generating track lines..."

timer = Timer()

parts = multiprocessing.cpu_count() - 1

prog = Progress(1 + parts + 2) #split, each sub-process, merge, delete temp

tally = 0

tmp_wkspc_path = os.path.split(sOutWorkspace)[0]

#split the point data into separate dictionaries for each process.

split_dicts = SplitDataDictEqually(dataDict, parts, iTotPoints)

tally+=1

prog.Update(tally)

#start each separate process asynchronously, and wait for all to finish before moving on.

pool = multiprocessing.Pool(parts)

jobs = []

for i in range(1,parts+1):

tmp_fgdb = os.path.join(tmp_wkspc_path, "temp" + str(i) + ".gdb")

sTempFC = tmp_fgdb+"\\"+sOutName+str(i)

splitDict = split_dicts[i-1]

tally+=1

p_callback = partial(MultiProcess_Status,p=prog,t=tally)

if sLineMethod == "Maximum Time and Distance":

jobs.append(pool.apply_async(AddLines_Segmented, (splitDict, sTempFC,iMaxTimeMinutes,iMaxDistMiles,spatialRef_CURSOR), callback=p_callback))

elif sLineMethod == "Time Interval":

jobs.append(pool.apply_async(AddLines_Interval, (splitDict, sTempFC,iIntervalHours,spatialRef_CURSOR), callback=p_callback))

pool.close()

pool.join()

del split_dicts

#Merge the temporary output feature classes in separate FGDBs into one feature class in the selected output geodatabase

to_merge = []

for i in range(1,parts+1):

tmp_fgdb = os.path.join(tmp_wkspc_path, "temp" + str(i) + ".gdb")

to_merge.append(tmp_fgdb+"\\"+sOutName+str(i))

arcpy.Merge_management(to_merge, sOutWorkspace+"\\"+sOutName)

tally+=1

prog.Update(tally)

#Delete the temporary geodatabases

for i in range(1,parts+1):

tmp_fgdb = os.path.join(tmp_wkspc_path, "temp" + str(i) + ".gdb")

arcpy.Delete_management(tmp_fgdb)

tally+=1

prog.Update(tally)

iTotTracks = int(arcpy.GetCount_management(sOutWorkspace+"\\"+sOutName).getOutput(0))

print "\r Unique "+sID_Field+"'s= "+FormatCounts(len(dataDict))

print " Track lines created = "+FormatCounts(iTotTracks)

print " Track lines created in:" + str(timer.End())

del timer

def MultiProcess_Status(x, p, t):

"""Simple function to get the status of each separate track line building process."""

p.Update(t)

def AddLines_Segmented(dataDict,sOutputFC,iMaxTimeMinutes,iMaxDistMiles,spatialRef_CURSOR,progress=False):

"""Creates track polylines and saves them in the ouput feature class. Track lines are created using a maximum distance and

maximum time threshold between points."""

prog = None

tally = 0

if progress:

prog = Progress(len(dataDict))

cur = arcpy.da.InsertCursor(sOutputFC, ("SHAPE@",sID_Field,"TrackStartTime","TrackEndTime"))

pnt1 = arcpy.Point()

pnt2 = arcpy.Point()

lineArray = arcpy.Array()

for key in sorted(dataDict.iterkeys()):

tally+=1

iSegCount = 0

#list of date time dictionary keys

dtTimes = sorted(dataDict[key].keys())

#start and end date variables: start date set with the start of any new line. end date constantly updated with the second point

#so that whenever the line ends the last datetime will already be set.

dtSegStart = None

dtSegEnd = None

#skip MMSI if there is only one point provided

if len(dtTimes) > 1:

for i in range(0,len(dtTimes)-1):

pnt1.X = dataDict[key][dtTimes[i]][0]

pnt1.Y = dataDict[key][dtTimes[i]][1]

pnt2.X = dataDict[key][dtTimes[i+1]][0]

pnt2.Y = dataDict[key][dtTimes[i+1]][1]

distance = GetDistance(pnt1.X,pnt1.Y,pnt2.X,pnt2.Y)

timeDiff = GetTimeDifference(dtTimes[i],dtTimes[i+1])

if distance < iMaxDistMiles and timeDiff < iMaxTimeMinutes:

if iSegCount == 0:

dtSegStart = dtTimes[i]

lineArray.add(pnt1)

lineArray.add(pnt2)

else:

lineArray.add(pnt2)

dtSegEnd = dtTimes[i+1]

iSegCount+=1

#Break the line as the gap exceeds the thresholds

else:

if lineArray.count > 1:

polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR)

cur.insertRow((polyline,key,dtSegStart,dtSegEnd))

del polyline

lineArray.removeAll()

dtSegStart = None

dtSegEnd = None

##reset seg count,since line ended early due to thresholds

iSegCount = 0

#end line since end of this MMSI set

if lineArray.count > 1:

polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR)

cur.insertRow((polyline,key,dtSegStart,dtSegEnd))

del polyline

lineArray.removeAll()

dtSegStart = None

dtSegEnd = None

if progress:

prog.Update(tally)

del prog

del cur

def AddLines_Interval(dataDict,sOutputFC,iIntervalHours,spatialRef_CURSOR,progress=False):

"""Creates track polylines and saves them in the ouput feature class. Track lines are created based on the specified time interval."""

prog = None

tally = 0

if progress:

prog = Progress(len(dataDict))

cur = arcpy.da.InsertCursor(sOutputFC, ("SHAPE@",sID_Field,"TrackStartTime","TrackEndTime"))

pnt1 = arcpy.Point()

pnt2 = arcpy.Point()

lineArray = arcpy.Array()

#create time interval from input parameter hours

tdInterval = datetime.timedelta(hours=iIntervalHours)

for key in sorted(dataDict.iterkeys()):

tally+=1

iSegCount = 0

#list of date time dictionary keys

dtTimes = sorted(dataDict[key].keys())

#start and end date variables: start date set with the start of any new line. end date constatnly updated with the second point

#so that whenever the line ends the last datetime will already be set.

dtSegStart = None

dtSegEnd = None

#initialize time interval using the first point as the start, and the end based on the selected interval

#updated as the interval threshold is reached.

dtIntervalBegin = dtTimes[0]

dtIntervalEnd = dtIntervalBegin + tdInterval

#skip MMSI if there is only one point provided

if len(dtTimes) > 1:

for i in range(0,len(dtTimes)-1):

pnt1.X = dataDict[key][dtTimes[i]][0]

pnt1.Y = dataDict[key][dtTimes[i]][1]

pnt2.X = dataDict[key][dtTimes[i+1]][0]

pnt2.Y = dataDict[key][dtTimes[i+1]][1]

if dtTimes[i] >= dtIntervalBegin and dtTimes[i+1] <= dtIntervalEnd:

if iSegCount == 0:

dtSegStart = dtTimes[i]

lineArray.add(pnt1)

lineArray.add(pnt2)

else:

lineArray.add(pnt2)

dtSegEnd = dtTimes[i+1]

iSegCount+=1

#Break the line as the next point is outside the defined interval

else:

if lineArray.count > 1:

polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR)

cur.insertRow((polyline,key,dtSegStart,dtSegEnd))

del polyline

lineArray.removeAll()

dtSegStart = None

dtSegEnd = None

#update for next time interval

dtIntervalBegin = dtIntervalEnd

dtIntervalEnd = dtIntervalBegin + tdInterval

##reset seg count,since line ended early due to thresholds

iSegCount = 0

#end line since end of this MMSI set

if lineArray.count > 1:

polyline = arcpy.Polyline(lineArray,spatialRef_CURSOR)

cur.insertRow((polyline,key,dtSegStart,dtSegEnd))

del polyline

lineArray.removeAll()

dtSegStart = None

dtSegEnd = None

dtIntervalBegin = None

dtIntervalEnd = None

if progress:

prog.Update(tally)

del prog

del cur

def AddVoyageData_NEW(dataDict,sOutputFC):

"""Adds voyage attribution from Voyage table to the assocatiated track lines. Voyage data is associated with the track lines based on the VoyageID.

Since one MMSI may have multiple VoyageIDs, the matching VoyageIDs can be selected by using the most frequent or the first."""

print "\n Adding Voyage data to output..."

timer = Timer()

prog = Progress(int(arcpy.GetCount_management(sOutputFC).getOutput(0)) + int(arcpy.GetCount_management(voyageTable).getOutput(0)))

tally = 0

iTotTracks = 0

iVoyageDataAdded = 0

#read ALL the voyage data into a python dictionary

dicVoyageAttrs = {}

with arcpy.da.SearchCursor(voyageTable,("VoyageID","Destination","Cargo","Draught","ETA","StartTime","EndTime")) as vRows:

for vRow in vRows:

tally+=1

dicVoyageAttrs[vRow[0]] = [vRow[1],vRow[2],vRow[3],vRow[4],vRow[5],vRow[6]]

prog.Update(tally)

#update the tracklines based on voyage data

with arcpy.da.UpdateCursor(sOutputFC, ("MMSI","TrackStartTime","TrackEndTime","VoyageID","Destination","Cargo","Draught","ETA","StartTime","EndTime")) as rows:

for row in rows:

tally+=1

iTotTracks+=1

lMMSI = row[0]

dtStart = row[1]

dtEnd = row[2]

if voyageMethod == "Most Frequent":

#Go through data dict and create dict of voyage IDs

dicVoyage = {}

for dt in sorted(dataDict[lMMSI].iterkeys()):

if dt >= dtStart and dt <= dtEnd:

iVoyageID = dataDict[lMMSI][dt][2]

if iVoyageID in dicVoyage:

dicVoyage[iVoyageID]+=1

else:

dicVoyage[iVoyageID] = 1

#get most frequent Voyage of Trackline

max_count = 0

for v in dicVoyage.iterkeys():

if dicVoyage[v] > max_count:

max_count = dicVoyage[v]

iVoyageID = v

else:#voyage method = "First"

for dt in sorted(dataDict[lMMSI].iterkeys()):

if dt >= dtStart and dt <= dtEnd:

iVoyageID = dataDict[lMMSI][dt][2]

break

#update row with voyageID

row[3] = iVoyageID

if dicVoyageAttrs.has_key(iVoyageID):

row[4] = dicVoyageAttrs[iVoyageID][0]

row[5] = dicVoyageAttrs[iVoyageID][1]

row[6] = dicVoyageAttrs[iVoyageID][2]

row[7] = dicVoyageAttrs[iVoyageID][3]

row[8] = dicVoyageAttrs[iVoyageID][4]

row[9] = dicVoyageAttrs[iVoyageID][5]

iVoyageDataAdded+=1

rows.updateRow(row)

prog.Update(tally)

print "\r Added to "+FormatCounts(iVoyageDataAdded)+" of "+FormatCounts(iTotTracks)+" track lines"

print " Added Voyage Data in:" + str(timer.End())

del prog

del timer

def AddVesselData_NEW(sOutputFC):

"""Adds vessel attribution from Vessel table to the assocatiated track lines. Vessel data is associated with the track lines based on the MMSI."""

print "\n Adding Vessel data to output..."

timer = Timer()

prog = Progress(int(arcpy.GetCount_management(sOutputFC).getOutput(0)) + int(arcpy.GetCount_management(vesselTable).getOutput(0)))

tally = 0

iTotTracks = 0

iVesselDataAdded = 0

#read ALL the vessel data into a python dictionary

dictVesselData = {}

with arcpy.da.SearchCursor(vesselTable,("MMSI","IMO","CallSign","Name","VesselType","Length","Width","DimensionComponents")) as vRows:

for vRow in vRows:

tally+=1

dictVesselData[vRow[0]] = [vRow[1],vRow[2],vRow[3],vRow[4],vRow[5],vRow[6],vRow[7]]

prog.Update(tally)

#update the tracklines based on vessel data

with arcpy.da.UpdateCursor(sOutputFC, ("MMSI","IMO","CallSign","Name","VesselType","Length","Width","DimensionComponents")) as rows:

for row in rows:

tally+=1

iTotTracks+=1

lMMSI = row[0]

if dictVesselData.has_key(lMMSI):

row[1] = dictVesselData[lMMSI][0]

row[2] = dictVesselData[lMMSI][1]

row[3] = dictVesselData[lMMSI][2]

row[4] = dictVesselData[lMMSI][3]

row[5] = dictVesselData[lMMSI][4]

row[6] = dictVesselData[lMMSI][5]

row[7] = dictVesselData[lMMSI][6]

iVesselDataAdded+=1

rows.updateRow(row)

prog.Update(tally)

print "\r Added to "+FormatCounts(iVesselDataAdded)+" of "+FormatCounts(iTotTracks)+" track lines"

print " Added Vessel Data in:" + str(timer.End())

del prog

del timer

def SetupVoyageCodedDomain(sOutputFC):

"""Setup the Coded Domains for the Cargo field. Creates a temporary table of the Cargo domain from the source database,

and then imports it to the ouput track feature class."""

print "\n Assigning Cargo coded domain values..."

input_wkspace = os.path.split(sInputFile)[0]

arcpy.env.workspace = sOutWorkspace

desc = arcpy.Describe(sOutWorkspace)

domains = desc.domains

#Check if domain already exists

bCargoDomExists = False

for domain in domains:

if domain == "Cargo":

bCargoDomExists = True

if bCargoDomExists == False:

#Copy the domain from the input AIS database to the output database

arcpy.DomainToTable_management(input_wkspace,"Cargo","CargoDom","code","description")

arcpy.TableToDomain_management("CargoDom","code","description",sOutWorkspace,"Cargo")

#Delete the temporary Domain Table

arcpy.Delete_management("CargoDom")

#Assign domain to the field in the track feature class

try:

arcpy.AssignDomainToField_management(sOutName,"Cargo","Cargo")

except:

print " ERROR Assigning Domain - Workspace is read only."

print " Domain can be assigned to the 'Cargo' field manually."

print " See the 'Use Limitations' section in the help documentation for details."

def SetupVesselCodedDomain(sOutputFC):

"""Setup the Coded Domains for the Vessel Type field. Creates a temporary table of the VesselType domain from the source database,

and then imports it to the ouput track feature class."""

print "\n Assigning Vessel Type coded domain values..."

input_wkspace = os.path.split(sInputFile)[0]

arcpy.env.workspace = sOutWorkspace

desc = arcpy.Describe(sOutWorkspace)

domains = desc.domains

#Check if domain already exists

bCargoDomExists = False

for domain in domains:

if domain == "VesselType":

bCargoDomExists = True

if bCargoDomExists == False:

#Copy the domain from the input AIS database to the output database

arcpy.DomainToTable_management(input_wkspace,"VesselType","VesselTypeDom","code","description")

arcpy.TableToDomain_management("VesselTypeDom","code","description",sOutWorkspace,"VesselType")

#Delete the temporary Domain Table

arcpy.Delete_management("VesselTypeDom")

#Assign domain to the field in the track feature class

try:

arcpy.AssignDomainToField_management(sOutName,"VesselType","VesselType")

except:

print " ERROR Assigning Domain - Workspace is read only."

print " Domain can be assigned to the 'VesselType' field manually."

print " See the 'Use Limitations' section in the help documentation for details."

def CreateTracks():

"""Main function that calls all of the other functions. Determines which functions are run, based on the user selected parameters."""

#Create Dictionary of broadcast data Key = ID, Value = Dictionary of {DateTime:[x,y,voyageID]}

dataDict, iTotPoints = CreateDataDict_NEW()

sOutputFC = sOutWorkspace+"\\"+sOutName

#If there are more than 2 CPUs then use multiprocessing approach, otherwise only use 1

# leaving one CPU to run operating system.

if multiprocessing.cpu_count() <= 2:

#Create output feature class

CreateOuputDataset(0)

#build track lines

print "\n Generating track lines..."

if sLineMethod == "Maximum Time and Distance":

AddLines_Segmented(dataDict,sOutputFC,iMaxTimeMinutes,iMaxDistMiles,spatialRef_CURSOR,True)

elif sLineMethod == "Time Interval":

AddLines_Interval(dataDict,sOutputFC,iIntervalHours,spatialRef_CURSOR,True)

iTotTracks = int(arcpy.GetCount_management(sOutputFC).getOutput(0))

print "\r Unique "+sID_Field+"'s= "+FormatCounts(len(dataDict))

print " Track lines created = "+FormatCounts(iTotTracks)

else:

#Create multiple temproary output feature classes

CreateOuputDataset(multiprocessing.cpu_count() - 1)

#Use multiprocessing method to create track lines

MultiProcess_Main(dataDict, iTotPoints)

#Add Voyage info, if table provided:

if sAddVoyageData == "true":

AddVoyageData_NEW(dataDict,sOutputFC)

SetupVoyageCodedDomain(sOutputFC)

#Add Vessel info, if table provided:

if sAddVesselData == "true":

AddVesselData_NEW(sOutputFC)

SetupVesselCodedDomain(sOutputFC)

print "\n"

def CheckParameters():

"""This function performs the checks required to make sure all the input parameters were provided and accurate.

These are the same checks performed within the Toolbox script tool validation code."""

global sAddVoyageData, sAddVesselData, spatialRef_CURSOR

global iMaxTimeMinutes, iMaxDistMiles, iIntervalHours

print "\n Checking Input Parameters..."

bPass = True

Messages = []

bRealAISdb = False

desc = None

#input points checks

if arcpy.Exists(sInputFile) == False:

bPass = False

Messages.append("Input Points do not exist")

else:

desc = arcpy.Describe(sInputFile)

if desc.datasetType <> "FeatureClass":

bPass = False

Messages.append("Input is not a feature class")

else:

if desc.shapeType <> "Point":

bPass = False

Messages.append("Input is not a point feature class")

else:

prjFile = os.path.join(arcpy.GetInstallInfo()["InstallDir"],"Coordinate Systems/Geographic Coordinate Systems/World/WGS 1984.prj")

spatialRef_WGS84 = arcpy.SpatialReference(prjFile)

spatialRef_input = arcpy.Describe(sInputFile).spatialReference

if spatialRef_WGS84.name <> spatialRef_input.name:

spatialRef_CURSOR = spatialRef_WGS84

else:

spatialRef_CURSOR = None

fields = desc.fields

#ID Field checks

bFieldExists = False

for field in fields:

if field.name == sID_Field:

bFieldExists = True

if field.type not in ["SmallInteger","Integer"]:

bPass = False

Messages.append("Input ID field is not an integer field")

break

if bFieldExists == False:

bPass = False

Messages.append("Input ID field does not exist in input feature class")

#DateTime Field checks

bFieldExists = False

for field in fields:

if field.name == sDT_Field:

bFieldExists = True

if field.type <> "Date":

bPass = False

Messages.append("Input Date/Time field is not an date field")

break

if bFieldExists == False:

bPass = False

Messages.append("Input Date/Time field does not exist in input feature class")

#check if adding vessel/voyage data is appropriate (real AIS database)

bRealAISdb = CheckAISDatabase()

#Output Workspace checks

if arcpy.Exists(sOutWorkspace) == False:

bPass = False

Messages.append("Output workspace does not exist")

else:

desc = arcpy.Describe(sOutWorkspace)

if desc.workspaceType <> "LocalDatabase":

bPass = False

Messages.append("Output workspace is not a geodatabase")

#Check output name

if sOutName <> str(arcpy.ValidateTableName(sOutName,sOutWorkspace)):

bPass = False

Messages.append("Output trackline feature class name is invalid")

#Check Break trackline option

if sLineMethod not in ["Maximum Time and Distance","Time Interval"]:

bPass = False

Messages.append("Break trackline method parameter is not a valid option")

#Check Maximum Time in minutes

try:

if sLineMethod == "Maximum Time and Distance":

test = int(sMaxTimeMinutes)

if int(sMaxTimeMinutes) <= 0:

bPass = False

Messages.append("Maximum Time must be greater than zero when using the Maximum Time and Distance option")

else:

iMaxTimeMinutes = int(sMaxTimeMinutes)

except:

bPass = False

Messages.append("Maximum Time is not numeric")

#Check Maximum Distance

try:

if sLineMethod == "Maximum Time and Distance":

test = int(sMaxDistMiles)

if int(sMaxDistMiles) <= 0:

bPass = False

Messages.append("Maximum Distance must be greater than zero when using the Maximum Time and Distance option")

else:

iMaxDistMiles = int(sMaxDistMiles)

except:

bPass = False

Messages.append("Maximum Distance is not numeric")

#Check Time Interval

try:

if sLineMethod == "Time Interval":

test = int(sIntervalHours)

if int(sIntervalHours) <= 0:

bPass = False

Messages.append("Time Interval must be greater than zero when using the Time Interval option")

else:

iIntervalHours = int(sIntervalHours)

except:

bPass = False

Messages.append("Time Interval is not numeric")

#Check Include voyage option

if sAddVoyageData not in ["true","false"]:

bPass = False

Messages.append("Option to include Voyage data must be 'true' or 'false'")

else:

if not bRealAISdb and sAddVoyageData == "true":

print " Ignoring Voyage data since input data is not part of a complete AIS database"

sAddVoyageData = "false"

elif sAddVoyageData == "true":

#check voyage table

if arcpy.Exists(voyageTable) == False:

bPass = False

Messages.append("Voyage table does not exist")

else:

desc = arcpy.Describe(voyageTable)

if desc.datasetType <> "Table":

bPass = False

Messages.append("Voyage table parameter is not a table")

#Check Voyage Matching Option

if voyageMethod not in ["Most Frequent","First"]:

bPass = False

Messages.append("Voyage Matching Option parameter is not a valid option")

#Check Include vessel data option

if sAddVesselData not in ["true","false"]:

bPass = False

Messages.append("Option to include Vessel data must be 'true' or 'false'")

else:

if not bRealAISdb and sAddVesselData == "true":

print " Ignoring Vessel data since input data is not part of a complete AIS database"

sAddVesselData = "false"

elif sAddVesselData == "true":

#check vessel table

if arcpy.Exists(vesselTable) == False:

bPass = False

Messages.append("Vessel table does not exist")

else:

desc = arcpy.Describe(vesselTable)

if desc.datasetType <> "Table":

bPass = False

Messages.append("Vessel table parameter is not a table")

return bPass,Messages, bRealAISdb

def HelpText(Messages = []):

"""Returns the command line syntax and examples of usage, when there was a mistake."""

if len(Messages) > 0:

print "\n\n"

for message in Messages:

print " ERROR: " + message

if len(sys.argv) <> 1:

print "\n\n"

print "Trackbuilder Command Line Syntax:"

print """"""

print ""

print "Examples:"

print """"C:\Python27\ArcGISx6410.3\python.exe" "C:\AIS\TrackBuilderVersion2.1_x64.py" "C:\AIS\Zone19_2011_07.gdb\Broadcast" "MMSI" "BaseDateTime" "C:\AIS\Zone19_2011_07.gdb" "TrackLine" "Maximum Time and Distance" "60" "6" "0" "true" "C:\AIS\Zone19_2011_07.gdb\Voyage" "Most Frequent" "true" "C:\AIS\Zone19_2011_07.gdb\Vessel" """

print ""

print """"C:\Python27\ArcGISx6410.3\python.exe" "C:\AIS\TrackBuilderVersion2.1_x64.py" "C:\AIS\Zone19_2011_07.gdb\Broadcast" "MMSI" "BaseDateTime" "C:\AIS\Zone19_2011_07.gdb" "TrackLine" "Time Interval" "0" "0" "24" "false" "" "" "false" "" """

##########################################

##########################################

if __name__ == '__main__':

timer = Timer()

bRun = True

Errors = []

global bRealAISdb

if len(sys.argv) == 1:

bRun, Errors, bRealAISdb = CheckParameters()

else:

#only proceed if all parameters are present

if len(sys.argv) <> 15 and len(sys.argv) <> 10: #all or without voyage/vessel options

bRun = False

else:

sInputFile = sys.argv[1]

sID_Field = sys.argv[2]

sDT_Field = sys.argv[3]

sOutWorkspace = sys.argv[4]

sOutName = sys.argv[5]

sLineMethod = sys.argv[6]

sMaxTimeMinutes = sys.argv[7]

sMaxDistMiles = sys.argv[8]

sIntervalHours = sys.argv[9]

if len(sys.argv) == 15:

sAddVoyageData = sys.argv[10]

voyageTable = sys.argv[11]

voyageMethod = sys.argv[12]

sAddVesselData = sys.argv[13]

vesselTable = sys.argv[14]

else:

sAddVoyageData = "false"

voyageTable = ""

voyageMethod = ""

sAddVesselData = "false"

vesselTable = ""

bRun, Errors, bRealAISdb = CheckParameters()

#Run the process if all parameters were validated

if bRun:

#Updated to work with ArcGIS Desktop 10.1 or later

if GetVersion() in ["10.1","10.2","10.2.1","10.2.2","10.3","10.3.1","10.4","10.4.1","10.5"]:

if sys.exec_prefix.find("x64") > 0:

try:

arcpy.env.overwriteOutput = True

CreateTracks()

print "Track Lines complete!"

print " Process complete in " + str(datetime.timedelta(seconds=timer.End()))

del timer

except arcpy.ExecuteError:

for msg in range(0, arcpy.GetMessageCount()):

if arcpy.GetSeverity(msg) == 2:

print arcpy.GetMessage(msg)

except:

tb = sys.exc_info()[2]

for i in range(0,len(traceback.format_tb(tb))):

tbinfo = traceback.format_tb(tb)[i]

msg = " ERROR: \n " + tbinfo + " Error Info:\n " + str(sys.exc_info()[1])

print msg

else:

msg = " ERROR: The x64 version of the TrackBuilder can only be used with a 64bit version of Python"

print msg

else:

msg = " ERROR: TrackBuilder can only be used with ArcGIS 10.1 or later"

print msg

#show syntax and examples if parameters were missing or invalid.

else:

HelpText(Errors)

4.再利用核密度分析工具，选择合适的参数值，生成轨迹热力图

5.打开新生成的图像的属性表，打开符号系统，选择合适的色带

最终效果图如下

转载自：https://blog.csdn.net/qq_912917507/article/details/81099656