Project

General

Profile

Statistics
| Branch: | Tag: | Revision:

pycama / src / pycama / OutlinePlot.py @ 834:d989df597b80

History | View | Annotate | Download (13.7 KB)

1
#!/usr/bin/env python3
2
# -*- coding: utf-8 -*-
3

    
4
# Copyright 2016-2017 Maarten Sneep, KNMI
5
#
6
# Redistribution and use in source and binary forms, with or without
7
# modification, are permitted provided that the following conditions are met:
8
#
9
# 1. Redistributions of source code must retain the above copyright notice,
10
#    this list of conditions and the following disclaimer.
11
#
12
# 2. Redistributions in binary form must reproduce the above copyright notice,
13
#    this list of conditions and the following disclaimer in the documentation
14
#    and/or other materials provided with the distribution.
15
#
16
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
20
# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
23
# ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26

    
27
## @file OutlinePlot.py
28
#
29
# This file defines a subclass of the pycama.AnalysisAndPlot.AnalysisAndPlot class.
30
# This subclass extracts the granule outlines from the input data.
31
#
32
# @author Maarten Sneep
33

    
34
import math
35
import logging
36
from argparse import Namespace
37
import warnings
38
warnings.filterwarnings("ignore", category=FutureWarning)
39

    
40
import numpy as np
41
import netCDF4
42
import h5py
43

    
44
from .AnalysisAndPlot import AnalysisAndPlot
45

    
46
from .utilities import *
47

    
48
## Make the granule outlines.
49
#
50
#  Also includes netCDF write & read routines. Plotting should be possible from
51
#  a restored object.
52
#
53
class OutlinePlot(AnalysisAndPlot):
54
    ## Define a few variables to collect the latitude & longitudes and the irradiance used for processing.
55
    def setup(self, **kwargs):
56
        self.keys = []
57
        self.orbits = []
58
        self.irradiance = {}
59
        self.latitude = {}
60
        self.longitude = {}
61

    
62
    ## Extract the required information from the input data.
63
    #
64
    #  The actual data extraction is done in the pycama.Reader.Reader.read() method,
65
    #  by calling the pycama.File.File.outline() method. Here we copy the data gathered in these methods.
66
    #
67
    #  In addition the irradiance reference is copied and the orbit numbers are gathered.
68
    def process(self):
69
        self.keys = sorted(self.input_variables.outlines.keys())
70
        n_keys = len(self.keys)
71
        for i, key in enumerate(self.keys):
72
            if i > 0:
73
                self.progress(100*i/n_keys)
74
            self.latitude[key] = self.input_variables.outlines[key]['latitude']
75
            self.longitude[key] = self.input_variables.outlines[key]['longitude']
76
            self.irradiance[key] = self.input_variables.irradiance[key]
77

    
78
        self.orbits = sorted(list(set([int(key.split('_')[0]) for key in self.keys])))
79

    
80
    ## Write processed data to output netcdf file.
81
    #
82
    #  @param fname File to write to
83
    #  @param mode  Writing mode, defaults to append.
84
    #
85
    #  Write data (including extraction specific dimensions) to the group with
86
    #  the name given in the `storage_group_name` property ("`outlineplot_data`" for this class).
87
    #
88
    #  The granule metadata that is not stored by the pycama.EventPlot.EventPlot class is stored here.
89
    def dump(self, fname, mode='a'):
90
        compress={'compression':'gzip', 'compression_opts':3, 'shuffle':True, 'fletcher32':True}
91
        time_step = self.time_index_in_output
92
        with h5py.File(fname, 'a') as ref:
93
            # group
94
            try:
95
                grp = ref.create_group(self.storage_group_name)
96
                grp.attrs['comment'] = 'Outlines of granules'
97
            except:
98
                grp = ref[self.storage_group_name]
99

    
100
            for key in self.keys:
101
                try:
102
                    dim = grp.create_dataset('key_'+key, (len(self.latitude[key]), ), dtype=np.int32)
103
                    dim[:] = np.arange(len(self.latitude[key]), dtype=np.int32)
104
                except:
105
                    pass
106

    
107
                orbit = int(key.split('_')[0])
108

    
109
                try:
110
                    lat = grp.create_dataset('latitude_'+key,
111
                                             (len(self.latitude[key]),),
112
                                             dtype=np.float32,
113
                                             **compress)
114
                    lat.dims.create_scale(dim, 'key_'+key)
115
                    lat.dims[0].attach_scale(dim)
116
                    lat.attrs['standard_name'] = 'latitude'
117
                    lat.attrs['units'] = 'degrees_north'
118
                    lat.attrs['comment'] = 'latitude for key {0}'.format(key)
119
                    lat.attrs['time_index'] = time_step
120
                    lat.attrs['orbit'] = orbit
121
                    lat.attrs['irradiance'] = self.irradiance[key]
122
                except:
123
                    lat = grp['latitude_'+key]
124

    
125
                try:
126
                    lon = grp.create_dataset('longitude_'+key,
127
                                             (len(self.latitude[key]),),
128
                                             dtype=np.float32,
129
                                             **compress)
130
                    lon.dims.create_scale(dim, 'key_'+key)
131
                    lon.dims[0].attach_scale(dim)
132
                    lon.attrs['standard_name'] = 'longitude'
133
                    lon.attrs['units'] = 'degrees_east'
134
                    lon.attrs['comment'] = 'longitude for key {0}'.format(key)
135
                    lon.attrs['time_index'] = time_step
136
                    lon.attrs['orbit'] = orbit
137
                    lon.attrs['irradiance'] = self.irradiance[key]
138
                except:
139
                    lon = grp['longitude_'+key]
140

    
141
                lat[:] = self.latitude[key]
142
                lon[:] = self.longitude[key]
143

    
144
            for product in self.input_variables.event_counting.keys():
145
                try:
146
                    pgrp = grp.create_group(product)
147
                    pgrp.attrs['comment'] = 'Orbit metadata counting for {0}'.format(product)
148
                except:
149
                    pgrp = grp[product]
150

    
151
                for orbit_id in [k for k in self.input_variables.event_counting[product].keys() if not k.endswith('_events')]:
152
                    orbit_dict = self.input_variables.event_counting[product][orbit_id]
153

    
154
                    try:
155
                        ogrp = pgrp.create_group(orbit_id)
156
                        ogrp.attrs['comment'] = 'Orbit metadata counting for {0} in orbit {1:05d}'.format(product, self.input_variables.event_counting[product][orbit_id]['orbit'])
157
                    except:
158
                        ogrp = pgrp[orbit_id]
159

    
160
                    for k, v in self.input_variables.event_counting[product][orbit_id].items():
161
                        ogrp.attrs[k] = v
162

    
163
                    for k, v in self.input_variables.input_pointer[product][orbit_id].items():
164
                        ogrp.attrs[k] = ", ".join(v)
165

    
166
    ## Read processed data from the input file, for specified time index.
167
    #
168
    #  @param fname NetCDF file with input data.
169
    #  @param time_index Time slice to read.
170
    def ingest(self, fname, time_index, exclude=None):
171
        self.time_index_in_output = time_index
172
        with h5py.File(fname, 'r') as ref:
173
            try:
174
                grp = ref[self.storage_group_name]
175
            except:
176
                self.logger.warning("outline data not in '%s'.", os.path.basename(fname))
177
                return False
178

    
179
            self.keys = [v[4:] for v in grp.keys() if v.startswith('key_')]
180
            for key in self.keys:
181
                lat = grp['latitude_'+key]
182
                lon = grp['longitude_'+key]
183
                if lat.attrs['time_index'] != time_index:
184
                    self.keys = list(filter(lambda a: a != key, self.keys))
185
                    continue
186

    
187
                self.latitude[key] = lat[:]
188
                self.longitude[key] = lon[:]
189
                try:
190
                    self.irradiance[key] = lon.attrs['irradiance']
191
                    if isinstance(self.irradiance[key], bytes):
192
                        self.irradiance[key] = str(self.irradiance[key], 'utf-8')
193
                except KeyError:
194
                    self.logger.warning("Unknown irradiance for '%s'.", key)
195
                    self.irradiance[key] = "Unknown"
196

    
197
            self.orbits = sorted(list(set([int(key.split('_')[0]) for key in self.keys])))
198
            # fake a reader object when ingesting.
199
            self.input_variables = Namespace(event_counting={}, input_pointer={})
200

    
201
            for pname in grp.keys():
202
                pgrp = grp[pname]
203
                if isinstance(pgrp, h5py.Dataset):
204
                    continue
205
                if 'comment' not in pgrp.attrs:
206
                    continue
207
                self.input_variables.event_counting[pname] = {}
208
                self.input_variables.input_pointer[pname] = {}
209
                for oname in pgrp.keys():
210
                    ogrp = pgrp[oname]
211
                    if isinstance(ogrp, h5py.Dataset):
212
                        continue
213
                    orbit = int(oname.split('_')[1])
214
                    granule_start = oname.split('_')[2]
215
                    self.input_variables.event_counting[pname][oname] = {}
216
                    self.input_variables.input_pointer[pname][oname] = {}
217
                    self.input_variables.event_counting[pname][oname]['orbit'] = orbit
218
                    for k, v in ogrp.attrs.items():
219
                        if isinstance(v, str):
220
                            self.input_variables.input_pointer[pname][oname][k] = [v]
221
                        else:
222
                            self.input_variables.event_counting[pname][oname][k] = v
223
        return True
224

    
225
    ## Merge data into a combined dataset.
226
    #
227
    #  @param other The object to be added to self,
228
    #               also an instance of the pycama.OutlinePlot.OutlinePlot class.
229
    def __iadd__(self, other):
230
        for key in other.keys:
231
            self.latitude[key] = other.latitude[key]
232
            self.longitude[key] = other.longitude[key]
233
            self.irradiance[key] = other.irradiance[key]
234
        orbits = other.orbits
235
        orbits.extend(self.orbits)
236
        self.orbits = sorted(list(set(orbits)))
237

    
238
    ## Make a plot of a specified variable
239
    #
240
    #  @param varname The name of the variable to plot.
241
    #  @param figure  The matplotlib.figure.Figure instance to plot to.
242
    #  @param ax      The matplotlib.axes.Axes object to use. Default is None, in that case `ax = matplotlib.pyplot.gca()` shall be used.
243
    #
244
    #  @return `False` if no plot could be created, `True` otherwise.
245
    #
246
    #  All outlines are plotted to a single figure, using rotating colors.
247
    def plot(self, figure, ax=None, ax_idx=None, resolution='l', orbit=None, **kwargs):
248
        import matplotlib
249
        matplotlib.use('svg') # to allow running without X11
250
        matplotlib.rc('font', family='DejaVu Sans')
251
        import cartopy.crs as ccrs
252
        import matplotlib.pyplot as plt
253
        import matplotlib.path as mpath
254
        import matplotlib.cm as cm
255
        from matplotlib.colors import Normalize, LogNorm
256
        import matplotlib.ticker as ticker
257

    
258
        if orbit is not None:
259
            self.logger.info("Plotting granule outline for orbit {0}".format(orbit))
260
        else:
261
            self.logger.info("Plotting granule outlines")
262

    
263
        orbits = [int(key.split('_')[0]) for key in self.keys]
264

    
265
        if orbit is not None and (not is_sequence(orbit)):
266
            orbit = [orbit]
267

    
268
        if ax is not None:
269
            self.logger.error("For geolocation plots the axis object should not be supplied. Axis control can use the ax_idx parameter.")
270
            return False
271

    
272
        if ax_idx is None or len(ax_idx) != 3:
273
            ax_idx = (1,1,1)
274

    
275
        # background map
276
        pcar = ccrs.PlateCarree()
277
        m = ccrs.Mollweide(central_longitude=0)
278
        alpha=1.0
279
        ax = figure.add_subplot(ax_idx[0], ax_idx[1], ax_idx[2], projection=m)
280
        ax.coastlines(linewidth=0.5, color='k', resolution=('110m' if resolution == 'l' else '50m'), zorder=3)
281
        yticks = np.arange(-90, 91, 30)
282
        xticks = np.arange(-180, 181, 30)
283
        ax.gridlines(xlocs=xticks, ylocs=yticks, crs=pcar, linewidth=0.5, color='k', zorder=4)
284
        mplotlib_colors = self.color_list(len(self.keys))
285

    
286
        i = 0
287
        for o, key in zip(orbits, self.keys):
288
            if orbit is None or o in orbit:
289
                lon = self.longitude[key]
290
                lat = self.latitude[key]
291
                dlon = list(np.fabs(lon[1:] - lon[:-1]))
292
                dlon.insert(0, 0.0)
293
                dlon = np.asarray(dlon)
294
                idx = np.arange(0, dlon.shape[0], 1)
295
                split_ranges = list(idx[dlon > 100])
296
                split_ranges.insert(0, 0)
297
                split_ranges.append(-1)
298
                for a,b in zip(split_ranges[0:-1], split_ranges[1:]):
299
                    ax.plot(lon[a:b], lat[a:b], label=key,
300
                            color=mplotlib_colors[i%(len(mplotlib_colors))],
301
                            linestyle='solid', transform=pcar)
302
                i += 1
303
        return True
304

    
305
    ## Generate a list of random colors
306
    #
307
    #  The colors are limited in length to avoid colors that do not contrast with the white background.
308
    #
309
    #  @param n Lenth of the list.
310
    def color_list(self, n):
311
        rval = []
312
        for i in range(n):
313
            length = 1.0
314
            while length > 0.7:
315
                c = np.random.uniform(size=3)
316
                length = np.sqrt(np.sum(c**2))/np.sqrt(3)
317
            rval.append(tuple(c))
318
        return rval