1#!/usr/bin/env python3
2
3"""GeoPlot widget.
4
5Requires `matplotlib` toolkit `basemap`.
6"""
7
8import logging
9from datetime import datetime, timedelta
10from collections import OrderedDict
11import csv
12
13# always import matplotlib_agg before matplotlib
14from chart.common import matplotlib_agg # (unused import) pylint: disable=W0611
15import numpy as np
16from PIL import Image
17from matplotlib import figure
18from matplotlib import cm
19from matplotlib import pyplot
20from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
21from matplotlib.collections import PolyCollection
22
23try:
24 from mpl_toolkits import basemap
25
26except ImportError:
27 logging.warning('Failed to import basemap library')
28 basemap = None
29
30from chart.common.path import Path
31from chart.db.connection import db_connect
32from chart.db import ts
33from chart.project import settings
34from chart.events.db import find_events
35from chart.reports.widget import Widget
36from chart.common.util import timedelta_mul
37from chart.plots.geoloc import Geoloc
38from chart.plots.geoloc import CannotGeolocate
39from chart.common.exceptions import ConfigError
40from chart.common.texttime import texttime_to_datetime
41from chart.plotviewer.make_url import make_plot_url
42from chart.reports.widget import WidgetOptionChoice
43
44db_conn = db_connect('FDF_EVENTS')
45
46IMG_EXT = '.png'
47HTML_EXT = '.html'
48
49logger = logging.getLogger()
50
51
52def get_mask_coordinates(mask_filename):
53 """Return coordinates to display an overlay over the Geoplot.
54
55 Read the coordinates from a CSV file and return them as lat/long lists.
56 """
57 lats = []
58 longs = []
59 with settings.REPORT_TEMPLATE_DIR.joinpath(mask_filename).open(mode='r') as csvfile:
60 mask_reader = csv.reader(csvfile)
61 next(mask_reader, None) # ignore header
62 for row in mask_reader:
63 lats.append(row[0])
64 longs.append(row[1])
65
66 return lats, longs
67
68
69class GeoPlotWidget(Widget):
70 """Geolocated plot of a single parameter."""
71
72 name = 'geoplot'
73 image = 'widgets/geo.png'
74 thumbnail = 'widgets/geo_sm.png'
75
76 url = 'http://tctrac/projects/chart/wiki/GeoplotWidget'
77
78 options = OrderedDict([
79 ('filename', {'type': 'string',
80 'optional': True}),
81 ('absolute-start-time', {
82 'description': 'Use fixed value such as "launch" or an ISO8601 data as the '
83 'start of the graph',
84 'type': 'datetime',
85 'optional': True}),
86 ('datapoint', {'type': 'field',
87 'optional': True,
88 'description': 'If present, timeseries data is plotted onto the image'}),
89 ('event', {'description': ('If present, events of the given class are plotted onto '
90 'the image'),
91 'multiple': True,
92 'type': OrderedDict([
93 ('eventname', { # note, it should be "eventclass" but I don't want to break
94 # existing templates
95 'type': 'eventclass'}),
96 ('hide-running-edac', {
97 'type': 'boolean',
98 'default': False}),
99 ('colour', {
100 'type': 'colour',
101 'default': '#0000ff'})])}),
102 ('width', {'type': 'uint',
103 'default': 1000,
104 'help': 'Image width',
105 'unit': 'px'}),
106 ('height', {'type': 'uint',
107 'optional': True,
108 'help': 'Image height. If omitted the image will be scaled according '
109 'to the width',
110 'unit': 'px'}),
111 ('thumbnail-width', {'type': 'uint',
112 'default': 680,
113 'help': 'Thumbnail image width',
114 'unit': 'px'}),
115 ('thumbnail-height', {'type': 'uint',
116 'optional': True,
117 'help': 'Thumbnail image height. If omitted the image will be scaled '
118 'according to the thumbnail width',
119 'unit': 'px'}),
120 ('projection', {'type': 'string',
121 'default': 'cyl',
122 'choices': [
123 WidgetOptionChoice(name='cyl',
124 description='Plain cylindrical projection'),
125 WidgetOptionChoice(name='ortho',
126 description='3d projection from infinite distance'),
127 ]}),
128 # geos, nsper :: satellite projections
129 ('blue-marble', {'type': 'boolean',
130 'default': False,
131 'label': 'Blue Marble',
132 'description': ('Use NASA Blue Marble background image instead of '
133 'flat colours for land and sea')}),
134 ('sea-colour', {'type': 'colour',
135 'label': 'Sea colour',
136 'default': '#f0f0f0'}),
137 ('land-colour', {'type': 'colour',
138 'label': 'Land colour',
139 'default': '#e0e0e0'}),
140 ('latitude', {'type': 'float',
141 'default': 0.0,
142 'description': 'Latitude (y-axis) centre point'}),
143 ('longitude', {'type': 'float',
144 'default': 0.0,
145 'description': 'Longitude (x-axis) centre point'}),
146 ('condition', {'type': 'string',
147 'optional': True}),
148 ('phase', {'type': 'string',
149 'optional': True,
150 'choices': [WidgetOptionChoice(name='ascending'),
151 WidgetOptionChoice(name='descending')]}),
152 ('marker-size', {'type': 'uint',
153 'default': 6,
154 'description': 'Size of marker used for plotting points'}),
155 ('max-duration', {'type': 'duration',
156 'optional': True,
157 'description':
158 'Allow image to be clipped to a certain time duration'}),
159 ('colourbar', {'type': bool,
160 'default': True,
161 'description':
162 'Show a colour bar'}),
163 ('overlay', {'description': 'Display an overlay over the Geoplot',
164 'multiple': True,
165 'type': OrderedDict([
166 ('mask', {
167 'type': str}),
168 ('transparency', {
169 'type': 'float',
170 'default': 0.15}),
171 ('colour', {
172 'type': str,
173 })])}),
174 ('show-empty', {'type': bool,
175 'default': True,
176 'description':
177 'Display the plot even when no data is available'}),
178 ('title', {'type': str,
179 'optional': True,
180 'description': 'Specify the title, otherwise use default title test'}),
181 ])
182
183 document_options = OrderedDict([
184 ('sid', {'type': 'sid'}),
185 ('sensing_start', {'type': 'datetime'}),
186 ('sensing_stop', {'type': 'datetime'})])
187
188 def __init__(self):
189 super(GeoPlotWidget, self).__init__()
190 # cache retrieved events so we don't requery them for the thumbnail and full sized images
191 self.events = None
192 self.title = None
193
194 def pre_html(self, document):
195 """Set our title attribute so the TOC widget can include it."""
196 c = self.config
197 if 'title' in c:
198 self.title = c['title']
199
200 elif 'datapoint' in c:
201 self.title = '{name}: {desc}'.format(
202 name=c['datapoint'].name, desc=c['datapoint'].description)
203
204 elif 'event' in c:
205 self.title = ', '.join(e['eventname'].name for e in c['event']) + ' events'
206
207 else:
208 self.title = ''
209
210 document.figures.append(self.title)
211
212 def get_filename(self, filename, datapoint, event, document):
213 """Create a suitable filename for the full sized image.
214 If `filename` is given, use a uniquified and countered version of that.
215 Otherwise look for either first `datapoint` or first `event` for a name.
216 Otherwise assign a incrementing number.
217 Finally add a suffix counter if the name is not unique within the report so far.
218 """
219 if filename is None:
220 if datapoint is not None:
221 filename = 'GEO_{name}{ext}'.format(name=datapoint.name, ext=IMG_EXT)
222
223 elif event is not None:
224 filename = 'GEO_' + '_'.join(e['eventname'].name for e in event) + IMG_EXT
225
226 else:
227 filename = 'GEO_{count}{ext}'.format(count=document.figure_count, ext=IMG_EXT)
228
229 # convert it to a report specific filename
230 filename = document.theme.mod_filename(document, filename)
231
232 # check for duplicates
233 # logging.debug('pre dup check filename ' + str(filename))
234 # for f in document.aux_files:
235 # logging.debug(f)
236 if filename in document.aux_files:
237 # logging.debug('incr')
238 suffix = 1
239 filename_without_ext = filename.stem # remove extension to the complete filename
240
241 while '{filename}_{suffix}{ext}'.format(filename=filename_without_ext,
242 suffix=suffix,
243 ext=IMG_EXT) in document.aux_files:
244 suffix += 1
245
246 filename = Path('{filename}_{suffix}{ext}'.format(
247 filename=filename_without_ext, suffix=suffix, ext=IMG_EXT))
248
249 return filename
250
251 def html(self, document):
252 """Render the geoplot."""
253
254 c = self.config
255 dc = document.config
256 sid = dc['sid']
257
258 if c['projection'] == 'cyl':
259 aspect = 2
260
261 elif c['projection'] == 'ortho':
262 aspect = 1
263
264 else:
265 logging.warning('Unknown projection {p}'.format(p=c['projection']))
266 aspect = 1
267
268 if 'absolute-start-time' in c:
269 if isinstance(c['absolute-start-time'], datetime):
270 sensing_start = c['absolute-start-time']
271
272 else:
273 sensing_start = texttime_to_datetime(c['absolute-start-time'], sid=sid)
274
275 if dc['sensing_stop'] < sensing_start:
276 sensing_stop = sensing_start + timedelta(seconds=1)
277
278 else:
279 sensing_start = dc['sensing_start']
280
281 sensing_stop = dc['sensing_stop']
282 if 'max-duration' in c:
283 new_sensing_stop = sensing_start + c['max-duration']
284 logging.info('Clipping stop time from {old} to {new}'.format(
285 old=sensing_stop, new=new_sensing_stop))
286 sensing_stop = new_sensing_stop
287
288 # basemap adds 19 pixels all around the edge
289 BASEMAP_MARGIN = 19
290 width = c['width'] - BASEMAP_MARGIN
291 height = c.get('height')
292 COLOURBAR_WIDTH = 50
293 if height is None:
294 height = c['width'] // aspect - COLOURBAR_WIDTH * c['colourbar']
295
296 height -= BASEMAP_MARGIN
297
298 logging.info('Full size dimensions {w} x {h}'.format(w=width, h=height))
299
300 # pick a filename
301 filename = self.get_filename(c.get('filename'), c.get('datapoint'), c.get('event'),
302 document)
303
304 datapoint_info, width, height = self.render(
305 sid, sensing_start, sensing_stop, width, height, c, filename)
306
307 thumbnail_filename = Path('{base}{thumb}{ext}'.format(
308 base=filename.stem,
309 ext=IMG_EXT,
310 thumb=document.theme.thumbnail_suffix))
311
312 thumbnail_width = c['thumbnail-width']
313 if 'thumbnail-height' in c:
314 thumbnail_height = c['thumbnail-height']
315
316 else:
317 thumbnail_height = thumbnail_width // aspect - COLOURBAR_WIDTH * c['colourbar']
318
319 thumbnail_height -= BASEMAP_MARGIN
320 thumbnail_width -= BASEMAP_MARGIN
321
322 _, thumbnail_width, thumbnail_height = self.render(
323 sid,
324 sensing_start,
325 sensing_stop,
326 thumbnail_width,
327 thumbnail_height,
328 c,
329 thumbnail_filename)
330
331 if 'datapoint' in c:
332 url = make_plot_url(
333 datapoints=[{'field': c['datapoint']}],
334 sid=sid,
335 sensing_start=sensing_start,
336 sensing_stop=sensing_stop,
337 plot_type='GEO',
338 calibrated=False)
339
340 else:
341 url = None
342
343 # only display Geoplot if there are samples to show
344 if c['show-empty'] is False:
345 show_plot = False
346 for item in datapoint_info:
347 if item['timeseries-points'] > 0:
348 show_plot = True
349
350 else:
351 show_plot = True
352
353 if show_plot:
354 document.append_figure(title=self.title,
355 live_url=url,
356 filename=thumbnail_filename,
357 width=thumbnail_width,
358 height=thumbnail_height,
359 zoom_filename=filename,
360 zoom_width=width,
361 zoom_height=height,
362 summary_info=None,
363 datapoint_info=datapoint_info,
364 widget=self)
365
366 def render(self, sid, sensing_start, sensing_stop, width, height, c, filename):
367 """Render a single image (either the thumbnail or full sized image) for `sid` from
368 `sensing_start` to `sensing_stop`, at size `width` by `height`, to `filename`.
369 Reads various config option from `c` which should probably be passed as separated
370 parameters."""
371 # note that the figsize parameter appears to be ignored
372 # do not call pylab.figure as that creates special Figure objects
373 # which globally cache all their data values
374 DPI = 100
375 # this is really bad - it leaks memory and the oo calls would be better - but
376 # I can't find a way to get colorbar() to work with the oo api
377 fig = pyplot.figure(figsize=(width / 100.0, height / 100),
378 # frameon=False,
379 dpi=100)
380 # fig = figure.Figure(figsize=(width / DPI, height / DPI), dpi=DPI)
381 # avoid a memory "leak" where matplotlib keeps copies of all image data
382 # FigureCanvas(fig)
383 # ax =
384 fig.add_axes((0, 0, 1, 1))
385 ax = fig.gca()
386
387 m = basemap.Basemap(projection=c['projection'],
388 fix_aspect=False,
389 lat_0=c['latitude'],
390 lon_0=c['longitude'],
391 # llcrnrlat=-90,
392 # urcrnrlat=90,
393 # llcrnrlon=-180,
394 # urcrnrlon=180,
395 resolution='c')
396 # ,
397 # ax=ax) # 'l', 'i', 'h', 'f'
398 # width, height, lon_0, lat_0
399
400 if c['blue-marble']:
401 m.bluemarble(scale=0.5) # 5400x2700
402
403 else:
404 m.drawcoastlines() # color
405 # draw parallels and meridians.
406 m.drawparallels(np.arange(-90., 91., 30.))
407 m.drawmeridians(np.arange(-180., 181., 60.))
408 # m.drawmapboundary(fill_color=(0.85, 0.85, 0.85)) # sea color
409 # m.fillcontinents(color=(0.65, 0.65, 0.65), # land color,
410 # lake_color=(0.85, 0.85, 0.85))
411 m.drawmapboundary(fill_color=c['sea-colour']) # sea color
412 m.fillcontinents(color=c['land-colour'], # land color,
413 lake_color=c['land-colour'])
414
415 datapoint_info = []
416
417 # this is a timeseries plot
418 if 'datapoint' in c:
419 # cax,
420 acc = sensing_start
421 BLOCK = timedelta(days=1)
422 ts_count = 0
423 while acc <= sensing_stop:
424 ts_count += self.plot_ts(context=m,
425 sid=sid,
426 sensing_start=acc,
427 sensing_stop=min(sensing_stop, acc + BLOCK),
428 datapoint=c['datapoint'],
429 phase=c.get('phase'),
430 condition=c.get('condition'),
431 marker_radius=c['marker-size'])
432 acc += BLOCK
433
434 datapoint_info.append({'field': c['datapoint'],
435 'timeseries-points': ts_count})
436
437 # draw a colour bar
438 # im = ax.imshow()
439 if c['colourbar']:
440 try:
441 bar = m.colorbar(location='right', pad='5%')
442 # bar = m.colorbar(ax=ax, cax=cax, location='right', pad='5%')
443 bar.set_label('{name} ({unit})'.format(
444 name=c['datapoint'].name, unit=c['datapoint'].unit))
445 except AttributeError as e:
446 logging.error(e)
447 except FloatingPointError as e:
448 # this happens if the data is all the same value because the matplotlib
449 # colorbar doesn't handle it
450 logging.error(e)
451
452 # display an overlay over Geoplot
453 if 'overlay' in c:
454 for ovl in c['overlay']:
455 lats, longs = get_mask_coordinates(ovl['mask'])
456 x, y = m(longs, lats)
457
458 # display single point on map
459 if len(lats) == 1:
460 m.scatter(x[0],
461 y[0],
462 30,
463 color=ovl['colour'],
464 marker='o',
465 alpha=ovl['transparency'])
466
467 # overlay a region of interest on the map
468 else:
469 data = np.array((x, y)).T
470 segs = [data, ]
471 poly = PolyCollection(segs,
472 alpha=ovl['transparency'],
473 antialiaseds=(1, ),
474 zorder=10)
475
476 poly.set_linewidth(1)
477 poly.set_facecolors(ovl['colour'])
478 poly.set_edgecolors(ovl['colour'])
479 ax.add_collection(poly)
480
481 # this is an events plot
482 if 'event' in c:
483 if len(c['event']) > 1:
484 raise ConfigError('GeoPlot widget can only handle a single event classname')
485
486 self.plot_events(context=m,
487 event_class=c['event'][0]['eventname'],
488 sid=sid,
489 sensing_start=sensing_start,
490 sensing_stop=sensing_stop,
491 colour=c['event'][0].get('colour'),
492 marker_radius=c['marker-size'],
493 hide_running_edac=c['event'][0]['hide-running-edac'])
494
495 # pyplot.title("Equidistant Cylindrical Projection")
496 # setting pad to zero removes the white padding from the edge of the image
497 # but also clips the border
498 fig.savefig(str(filename), bbox_inches='tight') # , pad_inches=0.0)
499
500 # read the actual width and height back from file as basemap doesn't create at the
501 # requested size
502 image = Image.open(str(filename))
503 width = image.size[0]
504 height = image.size[1]
505
506 # logging.info('Make geoplot size {w} x {h}'.format(w=width, h=height))
507
508 return datapoint_info, width, height
509
510 def plot_ts(self,
511 context,
512 sid,
513 sensing_start,
514 sensing_stop,
515 datapoint,
516 phase,
517 condition,
518 marker_radius):
519 """Render geolocated timeseries data."""
520 logging.info('Plotting subpoints from {strt} to {stop}'.format(
521 strt=sensing_start, stop=sensing_stop))
522 bind_vars = {'start_time': sensing_start,
523 'stop_time': sensing_stop}
524
525 and_clauses = ['sensing_time>=:start_time AND sensing_time<:stop_time']
526
527 if phase is not None:
528 anxs = db_conn.query('SELECT start_time '
529 'FROM fdf_events '
530 'WHERE scid=:sid_name '
531 'AND name=\'ANX\' '
532 'AND start_time>=:start_time AND start_time<:stop_time '
533 'ORDER BY start_time',
534 sid_name=sid.name,
535 start_time=sensing_start - timedelta(hours=8),
536 stop_time=sensing_stop + timedelta(hours=4)).fetchall()
537
538 if phase == 'ascending':
539 phase_start = 0.75
540 phase_stop = 1.25
541
542 elif phase == 'descending':
543 phase_start = 0.25
544 phase_stop = 0.75
545
546 else:
547 raise ConfigError('Phase must be either ascending or descending')
548
549 region_clauses = []
550 for i in range(len(anxs) - 1):
551 dur = anxs[i + 1][0] - anxs[i][0]
552 region_clauses.append(
553 '(sensing_time>=:start_{i} AND sensing_time<:stop_{i})'.format(
554 i=i))
555 bind_vars['start_{i}'.format(i=i)] = anxs[i][0] + timedelta_mul(dur, phase_start)
556 bind_vars['stop_{i}'.format(i=i)] = anxs[i][0] + timedelta_mul(dur, phase_stop)
557
558 if len(region_clauses) > 0:
559 and_clauses.append('(' + ' OR '.join(region_clauses) + ')')
560
561 logging.debug('Retrieve points from {regions} regions'.format(
562 regions=len(region_clauses)))
563
564 geoloc = Geoloc(sid, sensing_start, sensing_stop)
565
566 if condition is not None:
567 and_clauses.append(condition)
568
569 # For speed we could use np.fromiter instead of fetchall here
570 rows = ts.select(table=datapoint.table,
571 sensing_start=sensing_start,
572 sensing_stop=sensing_stop,
573 fields=('SENSING_TIME', datapoint),
574 sid=sid,
575 where=and_clauses,
576 bindvars=bind_vars).fetchall()
577
578 # logging.debug('Fetched {cc} rows'.format(cc=len(rows)))
579
580 lats = np.zeros((len(rows),), dtype=float)
581 lons = np.zeros((len(rows),), dtype=float)
582 data = np.zeros((len(rows),), dtype=float)
583 cc = 0
584 for i in range(len(rows)):
585 try:
586 lats[cc], lons[cc] = geoloc.lat_lon(rows[i][0])
587 except CannotGeolocate:
588 continue
589
590 data[cc] = rows[i][1]
591 cc += 1
592
593 lats = lats[:cc]
594 lons = lons[:cc]
595 data = data[:cc]
596
597 x, y = context(lons, lats)
598
599 # logging.debug('Plot final data sizes lats {lats} lons {lons} data {data} x
600 # {x} y {y}'.format(
601 # lats=lats.shape, lons=lons.shape, data=data.shape, x=x.shape, y=y.shape))
602
603 if len(data) > 0:
604 data += 1
605 cax = context.scatter(x,
606 y,
607 s=marker_radius, # size
608 c=data, # color
609 cmap=pyplot.cm.jet,
610 # cmap=pyplot.cm.autumn,
611 # norm=
612 zorder=2, # 1 placed data behind the land
613 linewidths=0)
614
615 # logging.debug('plotted {count} points'.format(count=data.size))
616 return data.size
617 # return cax, data.size
618
619 def plot_events(self,
620 context,
621 sid,
622 sensing_start,
623 sensing_stop,
624 event_class,
625 marker_radius,
626 colour,
627 hide_running_edac):
628 """Render geolocated CHART events."""
629 lons = []
630 lats = []
631 logger.debug('Seeking events for {sid} from {strt} to {stop} of type {typ}'.format(
632 sid=sid, strt=sensing_start, stop=sensing_stop, typ=event_class.name))
633
634 if self.events is None:
635 self.events = list(
636 find_events(
637 sid=sid,
638 start_time=sensing_start,
639 stop_time=sensing_stop,
640 event_class=event_class,
641 ordering='start_time'))
642
643 for e in self.events:
644 # logger.debug('Found {e}'.format(e=e))
645 geo = e.geoloc
646
647 # If the option to hide running edacs is selected, simply ignore them.
648 if hide_running_edac is True:
649 if e.instance_properties['delta_change'] == timedelta(0, 16):
650 continue
651
652 if e.instance_properties['delta_time'] > timedelta(0, 17):
653 continue
654
655 if geo is not None:
656 lons.append(geo['start_lon'])
657 lats.append(geo['start_lat'])
658
659 x, y = context(lons, lats)
660
661 if len(x) > 0:
662 context.scatter(x,
663 y,
664 s=marker_radius,
665 # c=1,
666 cmap=pyplot.cm.jet,
667 zorder=2,
668 linewidths=0,
669 color=colour)
670
671 return len(x)