schain Commit - r1509:3817490cb0ad · Jicamarca Repository

1

import os

1

import os

2

import datetime

2

import datetime

3

import numpy

3

import numpy

4

from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter

4

from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter

5

6

from schainpy.model.graphics.jroplot_base import Plot, plt

6

from schainpy.model.graphics.jroplot_base import Plot, plt

7

from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot

7

from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot

8

from schainpy.utils import log

8

from schainpy.utils import log

9

10

import wradlib.georef as georef

10

import wradlib.georef as georef

11

12

EARTH_RADIUS = 6.3710e3

12

EARTH_RADIUS = 6.3710e3

13

14

15

def ll2xy(lat1, lon1, lat2, lon2):

15

def ll2xy(lat1, lon1, lat2, lon2):

16

17

p = 0.017453292519943295

17

p = 0.017453292519943295

18

a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \

18

a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \

19

numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2

19

numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2

20

r = 12742 * numpy.arcsin(numpy.sqrt(a))

20

r = 12742 * numpy.arcsin(numpy.sqrt(a))

21

theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)

21

theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)

22

* numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))

22

* numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))

23

theta = -theta + numpy.pi/2

23

theta = -theta + numpy.pi/2

24

return r*numpy.cos(theta), r*numpy.sin(theta)

24

return r*numpy.cos(theta), r*numpy.sin(theta)

25

26

27

def km2deg(km):

27

def km2deg(km):

28

'''

28

'''

29

Convert distance in km to degrees

29

Convert distance in km to degrees

30

'''

30

'''

31

32

return numpy.rad2deg(km/EARTH_RADIUS)

32

return numpy.rad2deg(km/EARTH_RADIUS)

33

34

35

36

class SpectralMomentsPlot(SpectraPlot):

36

class SpectralMomentsPlot(SpectraPlot):

37

'''

37

'''

38

Plot for Spectral Moments

38

Plot for Spectral Moments

39

'''

39

'''

40

CODE = 'spc_moments'

40

CODE = 'spc_moments'

41

# colormap = 'jet'

41

# colormap = 'jet'

42

# plot_type = 'pcolor'

42

# plot_type = 'pcolor'

43

44

class DobleGaussianPlot(SpectraPlot):

44

class DobleGaussianPlot(SpectraPlot):

45

'''

45

'''

46

Plot for Double Gaussian Plot

46

Plot for Double Gaussian Plot

47

'''

47

'''

48

CODE = 'gaussian_fit'

48

CODE = 'gaussian_fit'

49

# colormap = 'jet'

49

# colormap = 'jet'

50

# plot_type = 'pcolor'

50

# plot_type = 'pcolor'

51

52

class DoubleGaussianSpectraCutPlot(SpectraCutPlot):

52

class DoubleGaussianSpectraCutPlot(SpectraCutPlot):

53

'''

53

'''

54

Plot SpectraCut with Double Gaussian Fit

54

Plot SpectraCut with Double Gaussian Fit

55

'''

55

'''

56

CODE = 'cut_gaussian_fit'

56

CODE = 'cut_gaussian_fit'

57

58

class SnrPlot(RTIPlot):

58

class SnrPlot(RTIPlot):

59

'''

59

'''

60

Plot for SNR Data

60

Plot for SNR Data

61

'''

61

'''

62

63

CODE = 'snr'

63

CODE = 'snr'

64

colormap = 'jet'

64

colormap = 'jet'

65

66

def update(self, dataOut):

66

def update(self, dataOut):

67

68

data = {

68

data = {

69

'snr': 10*numpy.log10(dataOut.data_snr)

69

'snr': 10*numpy.log10(dataOut.data_snr)

70

}

70

}

71

72

return data, {}

72

return data, {}

73

74

class DopplerPlot(RTIPlot):

74

class DopplerPlot(RTIPlot):

75

'''

75

'''

76

Plot for DOPPLER Data (1st moment)

76

Plot for DOPPLER Data (1st moment)

77

'''

77

'''

78

79

CODE = 'dop'

79

CODE = 'dop'

80

colormap = 'jet'

80

colormap = 'jet'

81

82

def update(self, dataOut):

82

def update(self, dataOut):

83

84

data = {

84

data = {

85

'dop': 10*numpy.log10(dataOut.data_dop)

85

'dop': 10*numpy.log10(dataOut.data_dop)

86

}

86

}

87

88

return data, {}

88

return data, {}

89

90

class PowerPlot(RTIPlot):

90

class PowerPlot(RTIPlot):

91

'''

91

'''

92

Plot for Power Data (0 moment)

92

Plot for Power Data (0 moment)

93

'''

93

'''

94

95

CODE = 'pow'

95

CODE = 'pow'

96

colormap = 'jet'

96

colormap = 'jet'

97

98

def update(self, dataOut):

98

def update(self, dataOut):

99

data = {

99

data = {

100

'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)

100

'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)

101

}

101

}

102

return data, {}

102

return data, {}

103

104

class SpectralWidthPlot(RTIPlot):

104

class SpectralWidthPlot(RTIPlot):

105

'''

105

'''

106

Plot for Spectral Width Data (2nd moment)

106

Plot for Spectral Width Data (2nd moment)

107

'''

107

'''

108

109

CODE = 'width'

109

CODE = 'width'

110

colormap = 'jet'

110

colormap = 'jet'

111

112

def update(self, dataOut):

112

def update(self, dataOut):

113

114

data = {

114

data = {

115

'width': dataOut.data_width

115

'width': dataOut.data_width

116

}

116

}

117

118

return data, {}

118

return data, {}

119

120

class SkyMapPlot(Plot):

120

class SkyMapPlot(Plot):

121

'''

121

'''

122

Plot for meteors detection data

122

Plot for meteors detection data

123

'''

123

'''

124

125

CODE = 'param'

125

CODE = 'param'

126

127

def setup(self):

127

def setup(self):

128

129

self.ncols = 1

129

self.ncols = 1

130

self.nrows = 1

130

self.nrows = 1

131

self.width = 7.2

131

self.width = 7.2

132

self.height = 7.2

132

self.height = 7.2

133

self.nplots = 1

133

self.nplots = 1

134

self.xlabel = 'Zonal Zenith Angle (deg)'

134

self.xlabel = 'Zonal Zenith Angle (deg)'

135

self.ylabel = 'Meridional Zenith Angle (deg)'

135

self.ylabel = 'Meridional Zenith Angle (deg)'

136

self.polar = True

136

self.polar = True

137

self.ymin = -180

137

self.ymin = -180

138

self.ymax = 180

138

self.ymax = 180

139

self.colorbar = False

139

self.colorbar = False

140

141

def plot(self):

141

def plot(self):

142

143

arrayParameters = numpy.concatenate(self.data['param'])

143

arrayParameters = numpy.concatenate(self.data['param'])

144

error = arrayParameters[:, -1]

144

error = arrayParameters[:, -1]

145

indValid = numpy.where(error == 0)[0]

145

indValid = numpy.where(error == 0)[0]

146

finalMeteor = arrayParameters[indValid, :]

146

finalMeteor = arrayParameters[indValid, :]

147

finalAzimuth = finalMeteor[:, 3]

147

finalAzimuth = finalMeteor[:, 3]

148

finalZenith = finalMeteor[:, 4]

148

finalZenith = finalMeteor[:, 4]

149

150

x = finalAzimuth * numpy.pi / 180

150

x = finalAzimuth * numpy.pi / 180

151

y = finalZenith

151

y = finalZenith

152

153

ax = self.axes[0]

153

ax = self.axes[0]

154

155

if ax.firsttime:

155

if ax.firsttime:

156

ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]

156

ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]

157

else:

157

else:

158

ax.plot.set_data(x, y)

158

ax.plot.set_data(x, y)

159

160

dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')

160

dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')

161

dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')

161

dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')

162

title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,

162

title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,

163

dt2,

163

dt2,

164

len(x))

164

len(x))

165

self.titles[0] = title

165

self.titles[0] = title

166

167

168

class GenericRTIPlot(Plot):

168

class GenericRTIPlot(Plot):

169

'''

169

'''

170

Plot for data_xxxx object

170

Plot for data_xxxx object

171

'''

171

'''

172

173

CODE = 'param'

173

CODE = 'param'

174

colormap = 'viridis'

174

colormap = 'viridis'

175

plot_type = 'pcolorbuffer'

175

plot_type = 'pcolorbuffer'

176

177

def setup(self):

177

def setup(self):

178

self.xaxis = 'time'

178

self.xaxis = 'time'

179

self.ncols = 1

179

self.ncols = 1

180

self.nrows = self.data.shape('param')[0]

180

self.nrows = self.data.shape('param')[0]

181

self.nplots = self.nrows

181

self.nplots = self.nrows

182

self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})

182

self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})

183

184

if not self.xlabel:

184

if not self.xlabel:

185

self.xlabel = 'Time'

185

self.xlabel = 'Time'

186

187

self.ylabel = 'Range [km]'

187

self.ylabel = 'Range [km]'

188

if not self.titles:

188

if not self.titles:

189

self.titles = ['Param {}'.format(x) for x in range(self.nrows)]

189

self.titles = ['Param {}'.format(x) for x in range(self.nrows)]

190

191

def update(self, dataOut):

191

def update(self, dataOut):

192

193

data = {

193

data = {

194

'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)

194

'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)

195

}

195

}

196

197

meta = {}

197

meta = {}

198

199

return data, meta

199

return data, meta

200

201

def plot(self):

201

def plot(self):

202

# self.data.normalize_heights()

202

# self.data.normalize_heights()

203

self.x = self.data.times

203

self.x = self.data.times

204

self.y = self.data.yrange

204

self.y = self.data.yrange

205

self.z = self.data['param']

205

self.z = self.data['param']

206

self.z = 10*numpy.log10(self.z)

206

self.z = 10*numpy.log10(self.z)

207

self.z = numpy.ma.masked_invalid(self.z)

207

self.z = numpy.ma.masked_invalid(self.z)

208

209

if self.decimation is None:

209

if self.decimation is None:

210

x, y, z = self.fill_gaps(self.x, self.y, self.z)

210

x, y, z = self.fill_gaps(self.x, self.y, self.z)

211

else:

211

else:

212

x, y, z = self.fill_gaps(*self.decimate())

212

x, y, z = self.fill_gaps(*self.decimate())

213

214

for n, ax in enumerate(self.axes):

214

for n, ax in enumerate(self.axes):

215

216

self.zmax = self.zmax if self.zmax is not None else numpy.max(

216

self.zmax = self.zmax if self.zmax is not None else numpy.max(

217

self.z[n])

217

self.z[n])

218

self.zmin = self.zmin if self.zmin is not None else numpy.min(

218

self.zmin = self.zmin if self.zmin is not None else numpy.min(

219

self.z[n])

219

self.z[n])

220

221

if ax.firsttime:

221

if ax.firsttime:

222

if self.zlimits is not None:

222

if self.zlimits is not None:

223

self.zmin, self.zmax = self.zlimits[n]

223

self.zmin, self.zmax = self.zlimits[n]

224

225

ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],

225

ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],

226

vmin=self.zmin,

226

vmin=self.zmin,

227

vmax=self.zmax,

227

vmax=self.zmax,

228

cmap=self.cmaps[n]

228

cmap=self.cmaps[n]

229

)

229

)

230

else:

230

else:

231

if self.zlimits is not None:

231

if self.zlimits is not None:

232

self.zmin, self.zmax = self.zlimits[n]

232

self.zmin, self.zmax = self.zlimits[n]

233

ax.collections.remove(ax.collections[0])

233

ax.collections.remove(ax.collections[0])

234

ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],

234

ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],

235

vmin=self.zmin,

235

vmin=self.zmin,

236

vmax=self.zmax,

236

vmax=self.zmax,

237

cmap=self.cmaps[n]

237

cmap=self.cmaps[n]

238

)

238

)

239

240

241

class PolarMapPlot(Plot):

241

class PolarMapPlot(Plot):

242

'''

242

'''

243

Plot for weather radar

243

Plot for weather radar

244

'''

244

'''

245

246

CODE = 'param'

246

CODE = 'param'

247

colormap = 'seismic'

247

colormap = 'seismic'

248

249

def setup(self):

249

def setup(self):

250

self.ncols = 1

250

self.ncols = 1

251

self.nrows = 1

251

self.nrows = 1

252

self.width = 9

252

self.width = 9

253

self.height = 8

253

self.height = 8

254

self.mode = self.data.meta['mode']

254

self.mode = self.data.meta['mode']

255

if self.channels is not None:

255

if self.channels is not None:

256

self.nplots = len(self.channels)

256

self.nplots = len(self.channels)

257

self.nrows = len(self.channels)

257

self.nrows = len(self.channels)

258

else:

258

else:

259

self.nplots = self.data.shape(self.CODE)[0]

259

self.nplots = self.data.shape(self.CODE)[0]

260

self.nrows = self.nplots

260

self.nrows = self.nplots

261

self.channels = list(range(self.nplots))

261

self.channels = list(range(self.nplots))

262

if self.mode == 'E':

262

if self.mode == 'E':

263

self.xlabel = 'Longitude'

263

self.xlabel = 'Longitude'

264

self.ylabel = 'Latitude'

264

self.ylabel = 'Latitude'

265

else:

265

else:

266

self.xlabel = 'Range (km)'

266

self.xlabel = 'Range (km)'

267

self.ylabel = 'Height (km)'

267

self.ylabel = 'Height (km)'

268

self.bgcolor = 'white'

268

self.bgcolor = 'white'

269

self.cb_labels = self.data.meta['units']

269

self.cb_labels = self.data.meta['units']

270

self.lat = self.data.meta['latitude']

270

self.lat = self.data.meta['latitude']

271

self.lon = self.data.meta['longitude']

271

self.lon = self.data.meta['longitude']

272

self.xmin, self.xmax = float(

272

self.xmin, self.xmax = float(

273

km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)

273

km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)

274

self.ymin, self.ymax = float(

274

self.ymin, self.ymax = float(

275

km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)

275

km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)

276

# self.polar = True

276

# self.polar = True

277

278

def plot(self):

278

def plot(self):

279

280

for n, ax in enumerate(self.axes):

280

for n, ax in enumerate(self.axes):

281

data = self.data['param'][self.channels[n]]

281

data = self.data['param'][self.channels[n]]

282

283

zeniths = numpy.linspace(

283

zeniths = numpy.linspace(

284

0, self.data.meta['max_range'], data.shape[1])

284

0, self.data.meta['max_range'], data.shape[1])

285

if self.mode == 'E':

285

if self.mode == 'E':

286

azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2

286

azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2

287

r, theta = numpy.meshgrid(zeniths, azimuths)

287

r, theta = numpy.meshgrid(zeniths, azimuths)

288

x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(

288

x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(

289

theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))

289

theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))

290

x = km2deg(x) + self.lon

290

x = km2deg(x) + self.lon

291

y = km2deg(y) + self.lat

291

y = km2deg(y) + self.lat

292

else:

292

else:

293

azimuths = numpy.radians(self.data.yrange)

293

azimuths = numpy.radians(self.data.yrange)

294

r, theta = numpy.meshgrid(zeniths, azimuths)

294

r, theta = numpy.meshgrid(zeniths, azimuths)

295

x, y = r*numpy.cos(theta), r*numpy.sin(theta)

295

x, y = r*numpy.cos(theta), r*numpy.sin(theta)

296

self.y = zeniths

296

self.y = zeniths

297

298

if ax.firsttime:

298

if ax.firsttime:

299

if self.zlimits is not None:

299

if self.zlimits is not None:

300

self.zmin, self.zmax = self.zlimits[n]

300

self.zmin, self.zmax = self.zlimits[n]

301

ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),

301

ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),

302

x, y, numpy.ma.array(data, mask=numpy.isnan(data)),

302

x, y, numpy.ma.array(data, mask=numpy.isnan(data)),

303

vmin=self.zmin,

303

vmin=self.zmin,

304

vmax=self.zmax,

304

vmax=self.zmax,

305

cmap=self.cmaps[n])

305

cmap=self.cmaps[n])

306

else:

306

else:

307

if self.zlimits is not None:

307

if self.zlimits is not None:

308

self.zmin, self.zmax = self.zlimits[n]

308

self.zmin, self.zmax = self.zlimits[n]

309

ax.collections.remove(ax.collections[0])

309

ax.collections.remove(ax.collections[0])

310

ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),

310

ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),

311

x, y, numpy.ma.array(data, mask=numpy.isnan(data)),

311

x, y, numpy.ma.array(data, mask=numpy.isnan(data)),

312

vmin=self.zmin,

312

vmin=self.zmin,

313

vmax=self.zmax,

313

vmax=self.zmax,

314

cmap=self.cmaps[n])

314

cmap=self.cmaps[n])

315

316

if self.mode == 'A':

316

if self.mode == 'A':

317

continue

317

continue

318

319

# plot district names

319

# plot district names

320

f = open('/data/workspace/schain_scripts/distrito.csv')

320

f = open('/data/workspace/schain_scripts/distrito.csv')

321

for line in f:

321

for line in f:

322

label, lon, lat = [s.strip() for s in line.split(',') if s]

322

label, lon, lat = [s.strip() for s in line.split(',') if s]

323

lat = float(lat)

323

lat = float(lat)

324

lon = float(lon)

324

lon = float(lon)

325

# ax.plot(lon, lat, '.b', ms=2)

325

# ax.plot(lon, lat, '.b', ms=2)

326

ax.text(lon, lat, label.decode('utf8'), ha='center',

326

ax.text(lon, lat, label.decode('utf8'), ha='center',

327

va='bottom', size='8', color='black')

327

va='bottom', size='8', color='black')

328

329

# plot limites

329

# plot limites

330

limites = []

330

limites = []

331

tmp = []

331

tmp = []

332

for line in open('/data/workspace/schain_scripts/lima.csv'):

332

for line in open('/data/workspace/schain_scripts/lima.csv'):

333

if '#' in line:

333

if '#' in line:

334

if tmp:

334

if tmp:

335

limites.append(tmp)

335

limites.append(tmp)

336

tmp = []

336

tmp = []

337

continue

337

continue

338

values = line.strip().split(',')

338

values = line.strip().split(',')

339

tmp.append((float(values[0]), float(values[1])))

339

tmp.append((float(values[0]), float(values[1])))

340

for points in limites:

340

for points in limites:

341

ax.add_patch(

341

ax.add_patch(

342

Polygon(points, ec='k', fc='none', ls='--', lw=0.5))

342

Polygon(points, ec='k', fc='none', ls='--', lw=0.5))

343

344

# plot Cuencas

344

# plot Cuencas

345

for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):

345

for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):

346

f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))

346

f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))

347

values = [line.strip().split(',') for line in f]

347

values = [line.strip().split(',') for line in f]

348

points = [(float(s[0]), float(s[1])) for s in values]

348

points = [(float(s[0]), float(s[1])) for s in values]

349

ax.add_patch(Polygon(points, ec='b', fc='none'))

349

ax.add_patch(Polygon(points, ec='b', fc='none'))

350

351

# plot grid

351

# plot grid

352

for r in (15, 30, 45, 60):

352

for r in (15, 30, 45, 60):

353

ax.add_artist(plt.Circle((self.lon, self.lat),

353

ax.add_artist(plt.Circle((self.lon, self.lat),

354

km2deg(r), color='0.6', fill=False, lw=0.2))

354

km2deg(r), color='0.6', fill=False, lw=0.2))

355

ax.text(

355

ax.text(

356

self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),

356

self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),

357

self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),

357

self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),

358

'{}km'.format(r),

358

'{}km'.format(r),

359

ha='center', va='bottom', size='8', color='0.6', weight='heavy')

359

ha='center', va='bottom', size='8', color='0.6', weight='heavy')

360

361

if self.mode == 'E':

361

if self.mode == 'E':

362

title = 'El={}$^\circ$'.format(self.data.meta['elevation'])

362

title = 'El={}$^\circ$'.format(self.data.meta['elevation'])

363

label = 'E{:02d}'.format(int(self.data.meta['elevation']))

363

label = 'E{:02d}'.format(int(self.data.meta['elevation']))

364

else:

364

else:

365

title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])

365

title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])

366

label = 'A{:02d}'.format(int(self.data.meta['azimuth']))

366

label = 'A{:02d}'.format(int(self.data.meta['azimuth']))

367

368

self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]

368

self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]

369

self.titles = ['{} {}'.format(

369

self.titles = ['{} {}'.format(

370

self.data.parameters[x], title) for x in self.channels]

370

self.data.parameters[x], title) for x in self.channels]

371

372

class WeatherParamsPlot(Plot):

372

class WeatherParamsPlot(Plot):

373

#CODE = 'RHI'

373

#CODE = 'RHI'

374

#plot_name = 'RHI'

374

#plot_name = 'RHI'

375

plot_type = 'scattermap'

375

plot_type = 'scattermap'

376

buffering = False

376

buffering = False

377

378

def setup(self):

378

def setup(self):

379

380

self.ncols = 1

380

self.ncols = 1

381

self.nrows = 1

381

self.nrows = 1

382

self.nplots= 1

382

self.nplots= 1

383

self.ylabel= 'Range [km]'

383

self.ylabel= 'Range [km]'

384

self.xlabel= 'Range [km]'

384

self.xlabel= 'Range [km]'

385

self.polar = True

385

self.polar = True

386

self.grid = True

386

self.grid = True

387

if self.channels is not None:

387

if self.channels is not None:

388

self.nplots = len(self.channels)

388

self.nplots = len(self.channels)

389

self.nrows = len(self.channels)

389

self.nrows = len(self.channels)

390

else:

390

else:

391

self.nplots = self.data.shape(self.CODE)[0]

391

self.nplots = self.data.shape(self.CODE)[0]

392

self.nrows = self.nplots

392

self.nrows = self.nplots

393

self.channels = list(range(self.nplots))

393

self.channels = list(range(self.nplots))

394

395

self.colorbar=True

395

self.colorbar=True

396

self.width =8

396

self.width =8

397

self.height =8

397

self.height =8

398

self.ini =0

398

self.ini =0

399

self.len_azi =0

399

self.len_azi =0

400

self.buffer_ini = None

400

self.buffer_ini = None

401

self.buffer_ele = None

401

self.buffer_ele = None

402

self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})

402

self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})

403

self.flag =0

403

self.flag =0

404

self.indicador= 0

404

self.indicador= 0

405

self.last_data_ele = None

405

self.last_data_ele = None

406

self.val_mean = None

406

self.val_mean = None

407

408

def update(self, dataOut):

408

def update(self, dataOut):

409

410

data = {}

410

data = {}

411

meta = {}

411

meta = {}

412

if hasattr(dataOut, 'dataPP_POWER'):

412

if hasattr(dataOut, 'dataPP_POWER'):

413

factor = 1

413

factor = 1

414

if hasattr(dataOut, 'nFFTPoints'):

414

if hasattr(dataOut, 'nFFTPoints'):

415

factor = dataOut.normFactor

415

factor = dataOut.normFactor

416

417

mask = dataOut.data_snr<self.snr_threshold

417

mask = dataOut.data_snr<self.snr_threshold

418

419

if 'pow' in self.attr_data[0].lower():

419

if 'pow' in self.attr_data[0].lower():

420

# data['data'] = 10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor))

420

# data['data'] = 10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor))

421

tmp = numpy.ma.masked_array(10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor)), mask=mask)

421

tmp = numpy.ma.masked_array(10*numpy.log10(10.0*getattr(dataOut, self.attr_data[0])/(factor)), mask=mask)

422

else:

422

else:

423

tmp = numpy.ma.masked_array(getattr(dataOut, self.attr_data[0]), mask=mask)

423

tmp = numpy.ma.masked_array(getattr(dataOut, self.attr_data[0]), mask=mask)

424

# tmp = getattr(dataOut, self.attr_data[0])

424

# tmp = getattr(dataOut, self.attr_data[0])

425

426

r = dataOut.heightList

426

r = dataOut.heightList

427

delta_height = r[1]-r[0]

427

delta_height = r[1]-r[0]

428

valid = numpy.where(r>=0)[0]

428

valid = numpy.where(r>=0)[0]

429

data['r'] = numpy.arange(len(valid))*delta_height

429

data['r'] = numpy.arange(len(valid))*delta_height

430

431

try:

431

try:

432

data['data'] = tmp[self.channels[0]][:,valid]

432

data['data'] = tmp[self.channels[0]][:,valid]

433

except:

433

except:

434

data['data'] = tmp[0][:,valid]

434

data['data'] = tmp[0][:,valid]

435

436

if dataOut.mode_op == 'PPI':

436

if dataOut.mode_op == 'PPI':

437

self.CODE = 'PPI'

437

self.CODE = 'PPI'

438

self.title = self.CODE

438

self.title = self.CODE

439

elif dataOut.mode_op == 'RHI':

439

elif dataOut.mode_op == 'RHI':

440

self.CODE = 'RHI'

440

self.CODE = 'RHI'

441

self.title = self.CODE

441

self.title = self.CODE

442

443

data['azi'] = dataOut.data_azi

443

data['azi'] = dataOut.data_azi

444

data['ele'] = dataOut.data_ele

444

data['ele'] = dataOut.data_ele

445

data['mode_op'] = dataOut.mode_op

445

data['mode_op'] = dataOut.mode_op

446

var = data['data'].flatten()

446

var = data['data'].flatten()

447

r = numpy.tile(data['r'], data['data'].shape[0]).reshape(data['data'].shape)*1000

447

r = numpy.tile(data['r'], data['data'].shape[0]).reshape(data['data'].shape)*1000

448

lla = georef.spherical_to_proj(r, data['azi'], data['ele'], (-75.295893, -12.040436, 3379.2147))

448

lla = georef.spherical_to_proj(r, data['azi'], data['ele'], (-75.295893, -12.040436, 3379.2147))

449

meta['lat'] = lla[:,:,1].flatten()[var.mask==False]

449

meta['lat'] = lla[:,:,1].flatten()[var.mask==False]

450

meta['lon'] = lla[:,:,0].flatten()[var.mask==False]

450

meta['lon'] = lla[:,:,0].flatten()[var.mask==False]

451

data['var'] = numpy.array([var[var.mask==False]])

451

data['var'] = numpy.array([var[var.mask==False]])

452

453

return data, meta

453

return data, meta

454

455

def plot(self):

455

def plot(self):

456

data = self.data[-1]

456

data = self.data[-1]

457

z = data['data']

457

z = data['data']

458

r = data['r']

458

r = data['r']

459

self.titles = []

459

self.titles = []

460

461

self.ymax = self.ymax if self.ymax else numpy.nanmax(r)

461

self.ymax = self.ymax if self.ymax else numpy.nanmax(r)

462

self.ymin = self.ymin if self.ymin else numpy.nanmin(r)

462

self.ymin = self.ymin if self.ymin else numpy.nanmin(r)

463

self.zmax = self.zmax if self.zmax else numpy.nanmax(z)

463

self.zmax = self.zmax if self.zmax else numpy.nanmax(z)

464

self.zmin = self.zmin if self.zmin is not None else numpy.nanmin(z)

464

self.zmin = self.zmin if self.zmin is not None else numpy.nanmin(z)

465

466

if data['mode_op'] == 'RHI':

466

if data['mode_op'] == 'RHI':

467

try:

467

try:

468

if self.data['mode_op'][-2] == 'PPI':

468

if self.data['mode_op'][-2] == 'PPI':

469

self.ang_min = None

469

self.ang_min = None

470

self.ang_max = None

470

self.ang_max = None

471

except:

471

except:

472

pass

472

pass

473

self.ang_min = self.ang_min if self.ang_min else 0

473

self.ang_min = self.ang_min if self.ang_min else 0

474

self.ang_max = self.ang_max if self.ang_max else 90

474

self.ang_max = self.ang_max if self.ang_max else 90

475

r, theta = numpy.meshgrid(r, numpy.radians(data['ele']) )

475

r, theta = numpy.meshgrid(r, numpy.radians(data['ele']) )

476

elif data['mode_op'] == 'PPI':

476

elif data['mode_op'] == 'PPI':

477

try:

477

try:

478

if self.data['mode_op'][-2] == 'RHI':

478

if self.data['mode_op'][-2] == 'RHI':

479

self.ang_min = None

479

self.ang_min = None

480

self.ang_max = None

480

self.ang_max = None

481

except:

481

except:

482

pass

482

pass

483

self.ang_min = self.ang_min if self.ang_min else 0

483

self.ang_min = self.ang_min if self.ang_min else 0

484

self.ang_max = self.ang_max if self.ang_max else 360

484

self.ang_max = self.ang_max if self.ang_max else 360

485

r, theta = numpy.meshgrid(r, numpy.radians(data['azi']) )

485

r, theta = numpy.meshgrid(r, numpy.radians(data['azi']) )

486

487

self.clear_figures()

487

self.clear_figures()

488

489

for i,ax in enumerate(self.axes):

489

for i,ax in enumerate(self.axes):

490

491

if ax.firsttime:

491

if ax.firsttime:

492

ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))

492

ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))

493

ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)

493

ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)

494

if data['mode_op'] == 'PPI':

494

if data['mode_op'] == 'PPI':

495

ax.set_theta_direction(-1)

495

ax.set_theta_direction(-1)

496

ax.set_theta_offset(numpy.pi/2)

496

ax.set_theta_offset(numpy.pi/2)

497

498

else:

498

else:

499

ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))

499

ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))

500

ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)

500

ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)

501

if data['mode_op'] == 'PPI':

501

if data['mode_op'] == 'PPI':

502

ax.set_theta_direction(-1)

502

ax.set_theta_direction(-1)

503

ax.set_theta_offset(numpy.pi/2)

503

ax.set_theta_offset(numpy.pi/2)

504

505

ax.grid(True)

505

ax.grid(True)

506

if data['mode_op'] == 'RHI':

506

if data['mode_op'] == 'RHI':

507

len_aux = int(data['azi'].shape[0]/4)

507

len_aux = int(data['azi'].shape[0]/4)

508

mean = numpy.mean(data['azi'][len_aux:-len_aux])

508

mean = numpy.mean(data['azi'][len_aux:-len_aux])

509

if len(self.channels) !=1:

509

if len(self.channels) !=1:

510

self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]

510

self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]

511

else:

511

else:

512

self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]

512

self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]

513

elif data['mode_op'] == 'PPI':

513

elif data['mode_op'] == 'PPI':

514

len_aux = int(data['ele'].shape[0]/4)

514

len_aux = int(data['ele'].shape[0]/4)

515

mean = numpy.mean(data['ele'][len_aux:-len_aux])

515

mean = numpy.mean(data['ele'][len_aux:-len_aux])

516

if len(self.channels) !=1:

516

if len(self.channels) !=1:

517

self.titles = ['PPI {} at EL: {} CH {}'.format(self.self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]

517

self.titles = ['PPI {} at EL: {} CH {}'.format(self.self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]

518

else:

518

else:

519

self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]

519

self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages

             import os
             import datetime
             import numpy
             from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter
             from schainpy.model.graphics.jroplot_base import Plot, plt
             from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
             from schainpy.utils import log
             import wradlib.georef as georef
             EARTH_RADIUS = 6.3710e3
             def ll2xy(lat1, lon1, lat2, lon2):
                 p = 0.017453292519943295
                 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
                     numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
                 r = 12742 * numpy.arcsin(numpy.sqrt(a))
                 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
                                       * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
                 theta = -theta + numpy.pi/2
                 return r*numpy.cos(theta), r*numpy.sin(theta)
             def km2deg(km):
                 '''
                 Convert distance in km to degrees
                 '''
                 return numpy.rad2deg(km/EARTH_RADIUS)
             class SpectralMomentsPlot(SpectraPlot):
                 '''
                 Plot for Spectral Moments
                 '''
                 CODE = 'spc_moments'
                 # colormap = 'jet'
                 # plot_type = 'pcolor'
             class DobleGaussianPlot(SpectraPlot):
                 '''
                 Plot for Double Gaussian Plot
                 '''
                 CODE = 'gaussian_fit'
                 # colormap = 'jet'
                 # plot_type = 'pcolor'
             class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
                 '''
                 Plot SpectraCut with Double Gaussian Fit
                 '''
                 CODE = 'cut_gaussian_fit'
             class SnrPlot(RTIPlot):
                 '''
                 Plot for SNR Data
                 '''
                 CODE = 'snr'
                 colormap = 'jet'
                 def update(self, dataOut):
                     data = {
                         'snr': 10*numpy.log10(dataOut.data_snr)
                     }
                     return data, {}
             class DopplerPlot(RTIPlot):
                 '''
                 Plot for DOPPLER Data (1st moment)
                 '''
                 CODE = 'dop'
                 colormap = 'jet'
                 def update(self, dataOut):
                     data = {
                         'dop': 10*numpy.log10(dataOut.data_dop)
                     }
                     return data, {}
             class PowerPlot(RTIPlot):
                 '''
                 Plot for Power Data (0 moment)
                 '''
                 CODE = 'pow'
                 colormap = 'jet'
                 def update(self, dataOut):
                     data = {
                         'pow': 10*numpy.log10(dataOut.data_pow/dataOut.normFactor)
                     }
                     return data, {}
             class SpectralWidthPlot(RTIPlot):
                 '''
                 Plot for Spectral Width Data (2nd moment)
                 '''
                 CODE = 'width'
                 colormap = 'jet'
                 def update(self, dataOut):
                     data = {
                         'width': dataOut.data_width
                     }
                     return data, {}
             class SkyMapPlot(Plot):
                 '''
                 Plot for meteors detection data
                 '''
                 CODE = 'param'
                 def setup(self):
                     self.ncols = 1
                     self.nrows = 1
                     self.width = 7.2
                     self.height = 7.2
                     self.nplots = 1
                     self.xlabel = 'Zonal Zenith Angle (deg)'
                     self.ylabel = 'Meridional Zenith Angle (deg)'
                     self.polar = True
                     self.ymin = -180
                     self.ymax = 180
                     self.colorbar = False
                 def plot(self):
                     arrayParameters = numpy.concatenate(self.data['param'])
                     error = arrayParameters[:, -1]
                     indValid = numpy.where(error == 0)[0]
                     finalMeteor = arrayParameters[indValid, :]
                     finalAzimuth = finalMeteor[:, 3]
                     finalZenith = finalMeteor[:, 4]
                     x = finalAzimuth * numpy.pi / 180
                     y = finalZenith
                     ax = self.axes[0]
                     if ax.firsttime:
                         ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
                     else:
                         ax.plot.set_data(x, y)
                     dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
                     dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
                     title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
                                                                                                  dt2,
                                                                                                  len(x))
                     self.titles[0] = title
             class GenericRTIPlot(Plot):
                 '''
                 Plot for data_xxxx object
                 '''
                 CODE = 'param'
                 colormap = 'viridis'
                 plot_type = 'pcolorbuffer'
                 def setup(self):
                     self.xaxis = 'time'
                     self.ncols = 1
                     self.nrows = self.data.shape('param')[0]
                     self.nplots = self.nrows
                     self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
                     if not self.xlabel:
                         self.xlabel = 'Time'
                     self.ylabel = 'Range [km]'
                     if not self.titles:
                         self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
                 def update(self, dataOut):
                     data = {
                         'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
                     }
                     meta = {}
                     return data, meta
                 def plot(self):
                     # self.data.normalize_heights()
                     self.x = self.data.times
                     self.y = self.data.yrange
                     self.z = self.data['param']
                     self.z = 10*numpy.log10(self.z)
                     self.z = numpy.ma.masked_invalid(self.z)
                     if self.decimation is None:
                         x, y, z = self.fill_gaps(self.x, self.y, self.z)
                     else:
                         x, y, z = self.fill_gaps(*self.decimate())
                     for n, ax in enumerate(self.axes):
                         self.zmax = self.zmax if self.zmax is not None else numpy.max(
                             self.z[n])
                         self.zmin = self.zmin if self.zmin is not None else numpy.min(
                             self.z[n])
                         if ax.firsttime:
                             if self.zlimits is not None:
                                 self.zmin, self.zmax = self.zlimits[n]
                             ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
                                                    vmin=self.zmin,
                                                    vmax=self.zmax,
                                                    cmap=self.cmaps[n]
                                                    )
                         else:
                             if self.zlimits is not None:
                                 self.zmin, self.zmax = self.zlimits[n]
                             ax.collections.remove(ax.collections[0])
                             ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
                                                    vmin=self.zmin,
                                                    vmax=self.zmax,
                                                    cmap=self.cmaps[n]
                                                    )
             class PolarMapPlot(Plot):
                 '''
                 Plot for weather radar
                 '''
                 CODE = 'param'
                 colormap = 'seismic'
                 def setup(self):
                     self.ncols = 1
                     self.nrows = 1
                     self.width = 9
                     self.height = 8
                     self.mode = self.data.meta['mode']
                     if self.channels is not None:
                         self.nplots = len(self.channels)
                         self.nrows = len(self.channels)
                     else:
                         self.nplots = self.data.shape(self.CODE)[0]
                         self.nrows = self.nplots
                         self.channels = list(range(self.nplots))
                     if self.mode == 'E':
                         self.xlabel = 'Longitude'
                         self.ylabel = 'Latitude'
                     else:
                         self.xlabel = 'Range (km)'
                         self.ylabel = 'Height (km)'
                     self.bgcolor = 'white'
                     self.cb_labels = self.data.meta['units']
                     self.lat = self.data.meta['latitude']
                     self.lon = self.data.meta['longitude']
                     self.xmin, self.xmax = float(
                         km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
                     self.ymin, self.ymax = float(
                         km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
                     #  self.polar = True
                 def plot(self):
                     for n, ax in enumerate(self.axes):
                         data = self.data['param'][self.channels[n]]
                         zeniths = numpy.linspace(
 , self.data.meta['max_range'], data.shape[1])
                         if self.mode == 'E':
                             azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
                             r, theta = numpy.meshgrid(zeniths, azimuths)
                             x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
                                 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
                             x = km2deg(x) + self.lon
                             y = km2deg(y) + self.lat
                         else:
                             azimuths = numpy.radians(self.data.yrange)
                             r, theta = numpy.meshgrid(zeniths, azimuths)
                             x, y = r*numpy.cos(theta), r*numpy.sin(theta)
                         self.y = zeniths
                         if ax.firsttime:
                             if self.zlimits is not None:
                                 self.zmin, self.zmax = self.zlimits[n]
                             ax.plt = ax.pcolormesh(  # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
                                 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
                                 vmin=self.zmin,
                                 vmax=self.zmax,
                                 cmap=self.cmaps[n])
                         else:
                             if self.zlimits is not None:
                                 self.zmin, self.zmax = self.zlimits[n]
                             ax.collections.remove(ax.collections[0])
                             ax.plt = ax.pcolormesh(  # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
                                 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
                                 vmin=self.zmin,
                                 vmax=self.zmax,
                                 cmap=self.cmaps[n])
                         if self.mode == 'A':
                             continue
                         # plot district names
                         f = open('/data/workspace/schain_scripts/distrito.csv')
                         for line in f:
                             label, lon, lat = [s.strip() for s in line.split(',') if s]
                             lat = float(lat)
                             lon = float(lon)
                             # ax.plot(lon, lat, '.b', ms=2)
                             ax.text(lon, lat, label.decode('utf8'), ha='center',
                                     va='bottom', size='8', color='black')
                         # plot limites
                         limites = []
                         tmp = []
                         for line in open('/data/workspace/schain_scripts/lima.csv'):
                             if '#' in line:
                                 if tmp:
                                     limites.append(tmp)
                                 tmp = []
                                 continue
                             values = line.strip().split(',')
                             tmp.append((float(values[0]), float(values[1])))
                         for points in limites:
                             ax.add_patch(
                                 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
                         # plot Cuencas
                         for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
                             f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
                             values = [line.strip().split(',') for line in f]
                             points = [(float(s[0]), float(s[1])) for s in values]
                             ax.add_patch(Polygon(points, ec='b', fc='none'))
                         # plot grid
                         for r in (15, 30, 45, 60):
                             ax.add_artist(plt.Circle((self.lon, self.lat),
                                                      km2deg(r), color='0.6', fill=False, lw=0.2))
                             ax.text(
                                 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
                                 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
                                 '{}km'.format(r),
                                 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
                     if self.mode == 'E':
                         title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
                         label = 'E{:02d}'.format(int(self.data.meta['elevation']))
                     else:
                         title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
                         label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
                     self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
                     self.titles = ['{} {}'.format(
                         self.data.parameters[x], title) for x in self.channels]
             class WeatherParamsPlot(Plot):
                 #CODE = 'RHI'
                 #plot_name = 'RHI'
                 plot_type = 'scattermap'
                 buffering = False
                 def setup(self):
                     self.ncols = 1
                     self.nrows = 1
                     self.nplots= 1
                     self.ylabel= 'Range [km]'
                     self.xlabel= 'Range [km]'
                     self.polar = True
                     self.grid = True
                     if self.channels is not None:
                         self.nplots = len(self.channels)
                         self.nrows = len(self.channels)
                     else:
                         self.nplots = self.data.shape(self.CODE)[0]
                         self.nrows = self.nplots
                         self.channels = list(range(self.nplots))
                     self.colorbar=True
                     self.width   =8
                     self.height  =8
                     self.ini     =0
                     self.len_azi =0
                     self.buffer_ini  = None
                     self.buffer_ele   = None
                     self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
                     self.flag    =0
                     self.indicador= 0
                     self.last_data_ele = None
                     self.val_mean      = None
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     if hasattr(dataOut, 'dataPP_POWER'):
                         factor = 1
                     if hasattr(dataOut, 'nFFTPoints'):
                         factor = dataOut.normFactor
                     mask = dataOut.data_snr<self.snr_threshold
                     if 'pow' in self.attr_data[0].lower():
                         # data['data'] = 10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor))
-                        tmp = numpy.ma.masked_array(10*numpy.log10(getattr(dataOut, self.attr_data[0])/(factor)), mask=mask)
+                        tmp = numpy.ma.masked_array(10*numpy.log10(10.0*getattr(dataOut, self.attr_data[0])/(factor)), mask=mask)
                     else:
                         tmp = numpy.ma.masked_array(getattr(dataOut, self.attr_data[0]), mask=mask)
                         # tmp = getattr(dataOut, self.attr_data[0])
                     r = dataOut.heightList
                     delta_height = r[1]-r[0]
                     valid = numpy.where(r>=0)[0]
                     data['r'] = numpy.arange(len(valid))*delta_height
                     try:
                         data['data'] = tmp[self.channels[0]][:,valid]
                     except:
                         data['data'] = tmp[0][:,valid]
                     if dataOut.mode_op == 'PPI':
                         self.CODE = 'PPI'
                         self.title = self.CODE
                     elif dataOut.mode_op == 'RHI':
                         self.CODE = 'RHI'
                         self.title = self.CODE
                     data['azi'] = dataOut.data_azi
                     data['ele'] = dataOut.data_ele
                     data['mode_op'] = dataOut.mode_op
                     var = data['data'].flatten()
                     r = numpy.tile(data['r'], data['data'].shape[0]).reshape(data['data'].shape)*1000
                     lla = georef.spherical_to_proj(r, data['azi'], data['ele'], (-75.295893, -12.040436, 3379.2147))
                     meta['lat'] = lla[:,:,1].flatten()[var.mask==False]
                     meta['lon'] = lla[:,:,0].flatten()[var.mask==False]
                     data['var'] = numpy.array([var[var.mask==False]])
                     return data, meta
                 def plot(self):
                     data = self.data[-1]
                     z = data['data']
                     r = data['r']
                     self.titles = []
                     self.ymax = self.ymax if self.ymax else numpy.nanmax(r)
                     self.ymin = self.ymin if self.ymin else numpy.nanmin(r)
                     self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
                     self.zmin = self.zmin if self.zmin is not None else numpy.nanmin(z)
                     if data['mode_op'] == 'RHI':
                         try:
                             if self.data['mode_op'][-2] == 'PPI':
                                 self.ang_min = None
                                 self.ang_max = None
                         except:
                             pass
                         self.ang_min = self.ang_min if self.ang_min else 0
                         self.ang_max = self.ang_max if self.ang_max else 90
                         r, theta = numpy.meshgrid(r, numpy.radians(data['ele']) )
                     elif data['mode_op'] == 'PPI':
                         try:
                             if self.data['mode_op'][-2] == 'RHI':
                                 self.ang_min = None
                                 self.ang_max = None
                         except:
                             pass
                         self.ang_min = self.ang_min if self.ang_min else 0
                         self.ang_max = self.ang_max if self.ang_max else 360
                         r, theta = numpy.meshgrid(r, numpy.radians(data['azi']) )
                     self.clear_figures()
                     for i,ax in enumerate(self.axes):
                         if ax.firsttime:
                             ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
                             ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
                             if data['mode_op'] == 'PPI':
                                 ax.set_theta_direction(-1)
                                 ax.set_theta_offset(numpy.pi/2)
                         else:
                             ax.set_xlim(numpy.radians(self.ang_min),numpy.radians(self.ang_max))
                             ax.plt = ax.pcolormesh(theta, r, z, cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
                             if data['mode_op'] == 'PPI':
                                 ax.set_theta_direction(-1)
                                 ax.set_theta_offset(numpy.pi/2)
                         ax.grid(True)
                         if data['mode_op'] == 'RHI':
                             len_aux = int(data['azi'].shape[0]/4)
                             mean = numpy.mean(data['azi'][len_aux:-len_aux])
                             if len(self.channels) !=1:
                                 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]
                             else:
                                 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
                         elif data['mode_op'] == 'PPI':
                             len_aux = int(data['ele'].shape[0]/4)
                             mean = numpy.mean(data['ele'][len_aux:-len_aux])
                             if len(self.channels) !=1:
                                 self.titles = ['PPI {} at EL: {} CH {}'.format(self.self.labels[x], str(round(mean,1)), x) for x in range(self.nrows)]
                             else:
                                 self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]