View Javadoc
1   package au.gov.amsa.geo.adhoc;
2   
3   import java.io.BufferedOutputStream;
4   import java.io.File;
5   import java.io.FileNotFoundException;
6   import java.io.FileOutputStream;
7   import java.io.IOException;
8   import java.io.PrintStream;
9   import java.text.DecimalFormat;
10  import java.text.ParseException;
11  import java.text.SimpleDateFormat;
12  import java.util.Arrays;
13  import java.util.Date;
14  import java.util.List;
15  import java.util.TimeZone;
16  import java.util.concurrent.TimeUnit;
17  
18  import org.apache.log4j.Logger;
19  
20  import com.github.davidmoten.grumpy.core.Position;
21  
22  import au.gov.amsa.geo.Eez;
23  import au.gov.amsa.geo.ShapefileUtil;
24  import au.gov.amsa.geo.TimedPosition;
25  import au.gov.amsa.geo.distance.OperatorEffectiveSpeedChecker;
26  import au.gov.amsa.geo.model.SegmentOptions;
27  import au.gov.amsa.gt.Shapefile;
28  import au.gov.amsa.risky.format.BinaryFixes;
29  import au.gov.amsa.risky.format.BinaryFixesFormat;
30  import au.gov.amsa.risky.format.Downsample;
31  import au.gov.amsa.risky.format.Fix;
32  import au.gov.amsa.util.identity.MmsiValidator2;
33  import rx.Observable;
34  import rx.observables.GroupedObservable;
35  import rx.schedulers.Schedulers;
36  
37  public class DistanceTravelledInEezMain {
38  
39      private static final Logger log = Logger.getLogger(DistanceTravelledInEezMain.class);
40  
41      private enum Location {
42          IN, OUT, UNKNOWN;
43      }
44  
45      private static final class State {
46  
47          private static final long HOUR_MILLIS = TimeUnit.HOURS.toMillis(1);
48  
49          String date;
50  
51          int mmsi;
52          Fix fix;
53          Location location;
54          double distanceKm;
55  
56          double totalTimeMs;
57  
58          public double totalTimeHours() {
59              return totalTimeMs / HOUR_MILLIS;
60          }
61  
62          public String formattedDate() {
63              return DistanceTravelledInEezMain.formattedDate(date);
64          }
65          
66      }
67  
68      static String formattedDate(String date) {
69          SimpleDateFormat sdfIn = new SimpleDateFormat("yyyy-MM-dd");
70          sdfIn.setTimeZone(TimeZone.getTimeZone("UTC"));
71          try {
72              Date d = sdfIn.parse(date);
73              // Sabine's preferred date format for importing into SPSS
74              SimpleDateFormat sdfOut = new SimpleDateFormat("dd-MMM-yyyy");
75              sdfOut.setTimeZone(TimeZone.getTimeZone("UTC"));
76              return sdfOut.format(d);
77          } catch (ParseException e) {
78              throw new RuntimeException(e);
79          }
80      }
81      
82      public static void main(String[] args) throws FileNotFoundException, IOException {
83          System.out.println("running");
84  
85          File tracks = new File("/home/dxm/combinedSortedTracks");
86          long t = System.currentTimeMillis();
87          List<File> files = Arrays.asList(tracks.listFiles());
88          files.sort((a, b) -> a.getName().compareTo(b.getName()));
89          try (PrintStream out = new PrintStream(new BufferedOutputStream(new FileOutputStream("target/output.csv")))) {
90              out.println(Vessel.headings());
91              Observable //
92                      .from(files) //
93                      .filter(x -> x.getName().endsWith(".track.gz")) //
94                      .filter(x -> x.getName().startsWith("2019")) //
95                      .filter(x -> !x.getName().startsWith("2020-08")) //
96                      .flatMap(file -> {
97                          log.info(file);
98  
99                          // Note that the Shapefile objects are not thread-safe so we make new one for   
100                         // each file to enable parallel processing
101 
102                         // used for intersections with eez boundary
103                         Shapefile eezLine = Eez.loadEezLine();
104 
105                         // used for contains tests
106                         Shapefile eezPolygon = Eez.loadEezPolygon();
107                         return BinaryFixes //
108                                 .from(file, true, BinaryFixesFormat.WITH_MMSI) //
109                                 .subscribeOn(Schedulers.computation()) //
110                                 .filter(f -> MmsiValidator2.INSTANCE.isValid(f.mmsi())) //
111                                 .groupBy(fix -> fix.mmsi()) //
112                                 .flatMap(o -> calculateDistance(file, eezLine, eezPolygon, o));
113                     }, Runtime.getRuntime().availableProcessors()) //
114                     .filter(x -> x.state.fix != null) //
115                     .toBlocking() //
116                     .forEach(x -> out.println(x.line()));
117         }
118         System.out.println((System.currentTimeMillis() - t) + "ms");
119     }
120 
121     private static Observable<? extends Vessel> calculateDistance(File file, Shapefile eezLine, Shapefile eezPolygon,
122             GroupedObservable<Integer, Fix> o) {
123         return Observable.defer(() -> {
124             State state = new State();
125             state.date = file.getName().substring(0, file.getName().indexOf(".track.gz"));
126             state.mmsi = o.getKey();
127             state.location = Location.UNKNOWN;
128             return o //
129                     .compose(Downsample.minTimeStep(5, TimeUnit.MINUTES)) //
130                     .lift(new OperatorEffectiveSpeedChecker(
131                             SegmentOptions.builder().acceptAnyFixHours(480L).maxSpeedKnots(50).build()))
132                     .filter(check -> check.isOk()) //
133                     .map(check -> check.fix()) //
134                     .doOnNext(fix -> {
135                         //TODO unit test
136                         boolean inside = eezPolygon.contains(fix.lat(), fix.lon());
137                         Location location = inside ? Location.IN : Location.OUT;
138                         if (state.location != Location.UNKNOWN) {
139                             boolean crossed = state.location != location;
140                             if (crossed) {
141                                 TimedPosition point = ShapefileUtil.findRegionCrossingPoint(eezLine, state.fix, fix);
142                                 final double distance;
143                                 if (location == Location.IN) {
144                                     distance = distanceKm(fix.lat(), fix.lon(), point.lat, point.lon);
145                                 } else {
146                                     distance = distanceKm(state.fix.lat(), state.fix.lon(), point.lat, point.lon);
147                                 }
148                                 state.distanceKm += distance;
149                                 state.totalTimeMs += distance
150                                         / distanceKm(state.fix.lat(), state.fix.lon(), fix.lat(), fix.lon())
151                                         * (fix.time() - state.fix.time());
152                             } else if (location == Location.IN) {
153                                 state.distanceKm += distanceKm(state.fix.lat(), state.fix.lon(), fix.lat(), fix.lon());
154                                 state.totalTimeMs += fix.time() - state.fix.time();
155                             }
156                         }
157                         state.fix = fix;
158                         state.location = location;
159                     }) //
160                     .count() //
161                     .map(count -> new Vessel(count, state));
162         });
163     }
164 
165     static final class Vessel {
166         final long count;
167 
168         State state;
169 
170         double totalTimeHours() {
171             return state.totalTimeMs / TimeUnit.HOURS.toMillis(1);
172         }
173 
174         Vessel(long count, State state) {
175             this.count = count;
176             this.state = state;
177         }
178 
179         static String headings() {
180             return "dateUTC,mmsi,aisClass,numReports,distanceNmInEEZ,elapsedTimeHoursInEEZ";
181         }
182 
183         String line() {
184             DecimalFormat df = new DecimalFormat("0.000");
185             return String.format("%s,%s,%s,%s,%s,%s", //
186                     state.formattedDate(), //
187                     state.mmsi, //
188                     state.fix.aisClass(), //
189                     count, //
190                     df.format(state.distanceKm / 1.852), //
191                     df.format(state.totalTimeHours()));
192         }
193 
194     }
195 
196     private static double distanceKm(double lat, double lon, double lat2, double lon2) {
197         return Position.create(lat, lon).getDistanceToKm(Position.create(lat2, lon2));
198     }
199 
200 }