pycama / src / pycama / OutlinePlot.py @ 834:d989df597b80
History  View  Annotate  Download (13.7 KB)
1 
#!/usr/bin/env python3


2 
# * coding: utf8 *

3  
4 
# Copyright 20162017 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], 'utf8') 
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
