# 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