442 lines
13 KiB
Plaintext
442 lines
13 KiB
Plaintext
# Script for importing AIS data (maritime trajectories from the Brest area in France)
|
|
#
|
|
# Required: folder SECONDO with data
|
|
# file Areas
|
|
# SymbolicTrajectoryAlgebra must be activated
|
|
#
|
|
# Preparations: place file Areas into secondo/bin directory
|
|
# enter location of SECONDO folder in the next line
|
|
#
|
|
# running times shown from Apple MacBook with SSD disk
|
|
|
|
|
|
|
|
# define the directory containing the data
|
|
let dir = '/home/ralf/Schreibtisch/Maritime/SECONDO/'
|
|
|
|
|
|
# import raw csv data
|
|
# the constant relation used in this command
|
|
# corresponds to the schema in the csv file
|
|
let Nari =
|
|
[const rel(tuple([
|
|
Mmsi: int,
|
|
Status: int,
|
|
Turn: real,
|
|
Speed: real,
|
|
Course: real,
|
|
Heading: int,
|
|
Lon: real,
|
|
Lat: real,
|
|
Time: longint
|
|
]))
|
|
value ()]
|
|
csvimport[dir + 'AIS Data/nari_dynamic.csv', 1, ""] consume
|
|
|
|
# about 6 minutes
|
|
# result size 19035630
|
|
|
|
# The ais data use as a time stamp an integer value
|
|
# meaning seconds to the start time used in Unix.
|
|
# For using in secondo, we have to convert this format
|
|
# into a datetime
|
|
|
|
let UnixStart = [const instant value "1970-1-1"]
|
|
|
|
|
|
# We convert the raw trajectory data into moving points.
|
|
# Here, we just use the positional data.
|
|
# For each Mmsi, a single trip is created from the
|
|
# sequence of (time, position) pairs. If there is
|
|
# a gap of 380 seconds between two reported positions,
|
|
# the trip contains also a gap.
|
|
|
|
let NariD = Nari feed
|
|
sortby[Mmsi, Time]
|
|
projectextend[Mmsi, Status, Turn, Speed, Course, Heading
|
|
; Pos: makepoint(.Lon, .Lat),
|
|
Time: UnixStart + create_duration(.Time, "s")]
|
|
groupby[Mmsi; Trip: group feed
|
|
approximate[Time, Pos, [const duration value (0 380000)] ]
|
|
]
|
|
consume
|
|
|
|
# about 2 - 3 minutes
|
|
# result size 5055
|
|
|
|
# Define the position of the antenna that is connected
|
|
# the the receiver recording the ais data.
|
|
|
|
let Antenna = [const point value (-4.57091 48.35923)]
|
|
|
|
# divide into trips with separation 2 hours
|
|
|
|
let NariTrips = NariD feed
|
|
projectextendstream[Mmsi
|
|
; Trip: .Trip sim_trips[ [const duration value (0 7200000)]] ]
|
|
consume
|
|
|
|
# 1:25 min
|
|
# result size 94281
|
|
|
|
# Here, also trips are created from an object trajectory.
|
|
# A counter for the trip number is added
|
|
# for each object.
|
|
|
|
let NariTrips2 = NariD feed
|
|
loopsel[fun(t: TUPLE)
|
|
attr(t, Trip) sim_trips[ [const duration value (0 7200000)]]
|
|
namedtransformstream[Trip]
|
|
addcounter[TripNo, 1]
|
|
extend[Mmsi: attr(t, Mmsi)]
|
|
project[Mmsi, TripNo, Trip]
|
|
]
|
|
consume
|
|
|
|
# 1:28 min
|
|
# result size 94281
|
|
|
|
# We simplify each trip using a variant of Douglas Peucker
|
|
# using a maximum derivation of 30 meters to the original
|
|
# trajectory.
|
|
# For checking the 'compression' rate, we add the number of
|
|
# units of the original trajectory and the simplified
|
|
# trajectory.
|
|
|
|
let NariTrips3 = NariTrips2 feed extend[Old: no_components(.Trip)]
|
|
projectextend[Mmsi, TripNo, Old
|
|
; Trip: simplify(.Trip, 30.0, [const geoid value WGS1984]) ]
|
|
extend[New: no_components(.Trip)]
|
|
consume
|
|
|
|
# 2:00 min
|
|
|
|
# Convert the simplified trips into a unit representation.
|
|
|
|
let NariUnits = NariTrips3 feed
|
|
projectextendstream[Mmsi; UTrip: units(.Trip)]
|
|
addcounter[UnitId, 1]
|
|
consume
|
|
|
|
# Old from NariD
|
|
# 4:01 min
|
|
# 16907598
|
|
|
|
# new
|
|
# 8.6 seconds
|
|
# result size 615453
|
|
|
|
# create a btree for easy finsing all units belonging
|
|
# to a certain mmsi
|
|
|
|
let NariUnits_Mmsi_btree = NariUnits createbtree[Mmsi]
|
|
|
|
# old
|
|
# 1.06 min
|
|
|
|
# new
|
|
# 3.9 seconds
|
|
|
|
|
|
# Import some non-moving data from a csv file.
|
|
|
|
let NariS =
|
|
[const rel(tuple([
|
|
SourceMMSI: int,
|
|
IMO: int,
|
|
CallSign: text,
|
|
ShipName: text,
|
|
ShipType: text,
|
|
ToBow: int,
|
|
ToStern: int,
|
|
ToStarBoard: int,
|
|
ToPort: int,
|
|
ETA: text,
|
|
Draught: real,
|
|
Destination: text,
|
|
MotherShipMMSI: int,
|
|
Time: longint
|
|
]))
|
|
value ()]
|
|
csvimport[dir + 'AIS Data/nari_static.csv', 1, ""] consume
|
|
|
|
# 42 seconds
|
|
# result size 1078617
|
|
|
|
# How many ships?
|
|
# Some mmsi are present multiple times in the csv file, e.g.,
|
|
# with a different destination, thus, we have to remove
|
|
# duplicates.
|
|
|
|
query NariS feedproject[SourceMMSI] sort rdup count
|
|
|
|
# 11 seconds
|
|
# result 4842
|
|
|
|
# Extract the pure static ship data.
|
|
let Ships = NariS feedproject[SourceMMSI, IMO, CallSign, ShipName, ShipType] sort rdup consume
|
|
|
|
# Create an index for the ships.
|
|
let Ships_SourceMMSI_btree = Ships createbtree[SourceMMSI]
|
|
|
|
|
|
|
|
# compute near collision situations
|
|
|
|
|
|
# A little test for building bounding boxes extended by 500 meters
|
|
|
|
query NariUnits feed extend[BBox: bbox(.UTrip), Box: bbox2d(.UTrip)]
|
|
extend[Box2: .Box extendGeo[500]]
|
|
extend[Box3: box3d(.Box2, deftime(.UTrip))]
|
|
head[3] consume
|
|
|
|
# For a near collision, the bounding boxes must not overlap.
|
|
# To exploit an r-tree for this purpose, we extend one of the
|
|
# bouding boxes by 500 meters.
|
|
query
|
|
NariUnits feed extend[Box: bbox(.UTrip)] {u1}
|
|
NariUnits feed extend[Box:
|
|
box3d(bbox2d(.UTrip) extendGeo[500], deftime(.UTrip))] {u2}
|
|
itSpatialJoin[Box_u1, Box_u2]
|
|
count
|
|
|
|
# 1:41 min
|
|
# result 3145061
|
|
|
|
# Matching units must come from different ships.
|
|
|
|
query
|
|
NariUnits feed extend[Box: bbox(.UTrip)] {u1}
|
|
NariUnits feed extend[Box:
|
|
box3d(bbox2d(.UTrip) extendGeo[500], deftime(.UTrip))] {u2}
|
|
itSpatialJoin[Box_u1, Box_u2]
|
|
filter[.Mmsi_u1 < .Mmsi_u2]
|
|
count
|
|
|
|
# 1:43 min
|
|
# result 737637
|
|
|
|
# Instead of only counting the candidates, we collect them
|
|
# into a relation here.
|
|
|
|
let Nearby =
|
|
NariUnits feed extend[Box: bbox(.UTrip)] {u1}
|
|
NariUnits feed extend[Box:
|
|
box3d(bbox2d(.UTrip) extendGeo[500], deftime(.UTrip))] {u2}
|
|
itSpatialJoin[Box_u1, Box_u2]
|
|
filter[.Mmsi_u1 < .Mmsi_u2]
|
|
consume
|
|
|
|
# 1:58 min
|
|
|
|
# Do the same thing as before.
|
|
# Here, a variant of the extendGeo operator is used instead of
|
|
# building the 3D-box manually.
|
|
|
|
let Nearby2 =
|
|
NariUnits feed extend[Box: bbox(.UTrip)] {u1}
|
|
NariUnits feed extend[Box: bbox(.UTrip) extendGeo[500]] {u2}
|
|
itSpatialJoin[Box_u1, Box_u2]
|
|
filter[.Mmsi_u1 < .Mmsi_u2]
|
|
consume
|
|
|
|
# Define the standard geod used in the data as an object.
|
|
let wgs1984 = create_geoid("WGS1984")
|
|
|
|
# Look at pairs of matching units. Select one by number.
|
|
# Reduce units to the overlap time.
|
|
# For testing, just take the one result.
|
|
|
|
query Nearby feed
|
|
extend[Speed1: val(final(speed(.UTrip_u1, wgs1984))) * 3.6 / 1.852,
|
|
Speed2: val(final(speed(.UTrip_u2, wgs1984))) * 3.6 / 1.852,
|
|
CommonTime: intersection(deftime(.UTrip_u1), deftime(.UTrip_u2))]
|
|
projectextend[
|
|
Speed1, Speed2; Mmsi1: .Mmsi_u1, Mmsi2: .Mmsi_u2,
|
|
UTrip1: (.UTrip_u1 atperiods .CommonTime) transformstream extract[Elem],
|
|
UTrip2: (.UTrip_u2 atperiods .CommonTime) transformstream extract[Elem]]
|
|
addcounter[No, 1]
|
|
filter[.No = 2]
|
|
head[1]
|
|
consume
|
|
|
|
# additionally consider direction and projection to space
|
|
|
|
query Nearby feed
|
|
extend[
|
|
Speed1: val(final(speed(.UTrip_u1, wgs1984))) * 3.6 / 1.852,
|
|
Speed2: val(final(speed(.UTrip_u2, wgs1984))) * 3.6 / 1.852,
|
|
CommonTime: intersection(deftime(.UTrip_u1), deftime(.UTrip_u2))]
|
|
projectextend[Speed1, Speed2; Mmsi1: .Mmsi_u1, Mmsi2: .Mmsi_u2,
|
|
UTrip1: (.UTrip_u1 atperiods .CommonTime) transformstream extract[Elem],
|
|
UTrip2: (.UTrip_u2 atperiods .CommonTime) transformstream extract[Elem]]
|
|
extend[Heading1: direction2heading(direction(val(initial(.UTrip1)),
|
|
val(final(.UTrip1))))]
|
|
extend[Heading2: direction2heading(direction(val(initial(.UTrip2)),
|
|
val(final(.UTrip2))))]
|
|
extend[Traj1: trajectory(.UTrip1)]
|
|
extend[Traj2: trajectory(.UTrip2)]
|
|
addcounter[No, 1]
|
|
filter[.No = 1]
|
|
head[1]
|
|
consume
|
|
|
|
# Importing environmental data
|
|
# adapt local directory
|
|
|
|
let AnchorageAreas =
|
|
dbimport2(dir+'/Anchorage Areas/Anchorage Areas.dbf')
|
|
shpimport2(dir + 'Anchorage Areas/Anchorage Areas.shp') namedtransformstream[Area]
|
|
obojoin
|
|
consume
|
|
|
|
let EuropeCoastline =
|
|
dbimport2(dir+'/European Coastline/Europe Coastline.dbf')
|
|
shpimport2(dir + '/European Coastline/Europe Coastline.shp' ) namedtransformstream[Coast]
|
|
obojoin
|
|
consume
|
|
|
|
let WorldSeas =
|
|
dbimport2(dir+'/IHO World Seas/World_Seas.dbf')
|
|
shpimport2(dir + '/IHO World Seas/World_Seas.shp' ) namedtransformstream[Sea]
|
|
obojoin
|
|
consume
|
|
|
|
let ports =
|
|
dbimport2(dir+'/Ports of Brittany/port.dbf')
|
|
shpimport2(dir + '/Ports of Brittany/port.shp' ) namedtransformstream[Port]
|
|
obojoin
|
|
consume
|
|
|
|
let tssQuessant =
|
|
dbimport2(dir+'/TSS Ouessant/TSS Ouessant.dbf')
|
|
shpimport2(dir + '/TSS Ouessant/TSS Ouessant.shp' ) namedtransformstream[Quessant]
|
|
obojoin
|
|
consume
|
|
|
|
|
|
# define another directory containing source data
|
|
let dir2 = dir + 'AIS Data/Status, Codes and Types/'
|
|
|
|
# Import some csv files.
|
|
# First define the schema used in the file as an empty relation
|
|
# and then actually import the data.
|
|
|
|
let MmsiCountryCodesT = [ const rel(tuple([Code : int, Country : text])) value ()]
|
|
|
|
let MmsiCountryCodes = MmsiCountryCodesT csvimport[dir2+'MMSI Country Codes.csv',1,""] consume
|
|
|
|
let NaviStatusT = [ const rel(tuple([Code : int, Status : text])) value ()]
|
|
|
|
let NaviStatus = NaviStatusT csvimport[dir2+'Navigational Status.csv',1,""] consume
|
|
|
|
let shipTypeT = [const rel(tuple([ID_shipType : int, Shiptype_min : int, ShipType_max : int,
|
|
Type_Name : text, Ais_type_summary : text]))
|
|
value ()]
|
|
|
|
let shipType = shipTypeT csvimport[dir2+'Ship Types List.csv',1,""] consume
|
|
|
|
let shipTypeDetailT = [const rel(tuple([Id_DetailedType : int, Detailed_Type : text,
|
|
Id_ShipType : int]))
|
|
value () ]
|
|
|
|
let shipTypeDetail = shipTypeDetailT csvimport[dir2+'Ship Types Detailled List.csv',1,""] consume
|
|
|
|
|
|
# Devide the coast line of Europe into single segments
|
|
|
|
let EuropeCoast = EuropeCoastline feed projectextendstream[; Seg: segments(.Coast)] consume
|
|
|
|
# 2:31 min
|
|
# result size 2847446
|
|
|
|
|
|
# Constructing a table of interesting areas
|
|
# constructed manually
|
|
#
|
|
# The following only works on the machine where these areas have been drawn manually.
|
|
|
|
# let Areas =
|
|
# Transrade feed namedtransformstream[Area] extend[Label: [const label value "Transrade"]]
|
|
# NavalAcademy feed namedtransformstream[Area] extend[Label: [const label value "NavalAcademy"]] concat
|
|
# Goulet feed namedtransformstream[Area] extend[Label: [const label value "Goulet"]] concat
|
|
# Camaret feed namedtransformstream[Area] extend[Label: [const label value "Camaret"]] concat
|
|
# BrestMilitary feed namedtransformstream[Area] extend[Label: [const label value "BrestMilitary"]] concat
|
|
# BrestSailing1 feed namedtransformstream[Area] extend[Label: [const label value "BrestSailing1"]] concat
|
|
# BrestSailing2 feed namedtransformstream[Area] extend[Label: [const label value "BrestSailing2"]] concat
|
|
# BrestCommercial feed namedtransformstream[Area] extend[Label: [const label value "BrestCommercial"]] concat
|
|
# Douarnenez feed namedtransformstream[Area] extend[Label: [const label value "Douarnenez"]] concat
|
|
# Morgat feed namedtransformstream[Area] extend[Label: [const label value "Morgat"]] concat
|
|
# Conquet feed namedtransformstream[Area] extend[Label: [const label value "Conquet"]] concat
|
|
# Out feed namedtransformstream[Area] extend[Label: [const label value "Out"]] concat
|
|
# consume
|
|
|
|
# save Areas to Areas
|
|
|
|
# on other machines, get the file Areas and put it into the bin directory. Then
|
|
|
|
restore Areas from Areas
|
|
|
|
# Extract a single object from the relation of areas.
|
|
let Transrade = Areas feed filter[.Label = "Transrade"] extract[Area]
|
|
|
|
# find all trips passing the Transrade Area
|
|
let R1 = NariTrips3 feed extend[Box: bbox2d(.Trip)]
|
|
Areas feed extend[Box2: bbox(.Area)]
|
|
itSpatialJoin[Box, Box2]
|
|
filter[.Trip passes .Area]
|
|
consume
|
|
|
|
|
|
# Collect the units of the trips that intersect an area of the Areas relation.
|
|
|
|
let R2 = R1 feed extendstream[Piece: components(.Trip at .Area)] consume
|
|
|
|
# Append the currently passed Area as a moving label.
|
|
let R3 = R2 feed extend[Time0: inst(initial(.Piece))]
|
|
extend[Unit: the_unit(.Label, getIntervals(deftime(.Piece)) transformstream extract[Elem])]
|
|
sortby[Mmsi, TripNo, Time0]
|
|
groupby[Mmsi, TripNo; Symbolic: group feed makemvalue[Unit]]
|
|
consume
|
|
|
|
|
|
|
|
# Name: matches
|
|
# Signature: {mlabel(s)|mplace(s)} x {pattern|text} -> bool
|
|
# Syntax: M matches P
|
|
# Meaning: Checks whether the trajectory M matches the pattern
|
|
# P.
|
|
# Example: query mlabel1 matches '(_ "Eving") * (_ "Innenstadt-Ost") +'
|
|
|
|
# R3 contains only the parts of a trip that are inside an area.
|
|
# Here, append the whole trip.
|
|
let NariTrips4 = NariTrips3 feed {n}
|
|
R3 feed
|
|
itHashJoin[Mmsi_n, Mmsi]
|
|
filter[.TripNo_n = .TripNo]
|
|
consume
|
|
|
|
|
|
# Here are some examples of pattern queries:
|
|
|
|
# find trips involving a subtrip from Transrade to NavalAcademy
|
|
|
|
# query NariTrips4 feed filter[.Symbolic matches '* (_ "Transrade") * (_ "NavalAcademy") *'] extend[Traj: trajectory3(.Trip_n)] consume
|
|
|
|
# find trips starting in Morgat
|
|
|
|
# query NariTrips4 feed filter[.Symbolic matches '(_ "Morgat") * '] extend[Traj: trajectory3(.Trip_n)] consume
|
|
|
|
# find trips only in Douarnenez
|
|
|
|
# query NariTrips4 feed filter[.Symbolic matches '(_ "Douarnenez")'] extend[Traj: trajectory3(.Trip_n)] consume
|
|
|
|
# find trips only in Douarnenez at night
|
|
|
|
# query NariTrips4 feed filter[.Symbolic matches '(night "Douarnenez")'] extend[Traj: trajectory3(.Trip_n)] consume
|
|
|
|
|
|
|
|
|