python 轨迹 车辆_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.&