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

You may also like...

退出移动版