ArcGIS+ArcPy制作船舶（车辆）轨迹热力图
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 """<path to x64 python.exe> <path to Trackbuilder .py script> <Path to input points> <ID field name> <DateTime field name> <Path to output workspace> <name of output trackline feature class> <Break Track line method option> <Maximum Time in minutes> <Maximum Distance in miles> <Time Interval in hours> <Include Voyage Attributes> <Voyage Table> <Voyage Matching Option> <Include Vessel Attributes> <Vessel Table>"""
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
