rad_count_dist with plotting utility

goesr/rad_count_dist.py

+#!/usr/bin/env python
+"""
+Compute distribution of unscaled pixel values in GOES-16 ABI L1B files
+Copyright 2017 University of Wisconsin Regents
+License: GPLv3 as of 2017-01-01
+Author: R.K.Garcia <rayg@ssec.wisc.edu>
+Recommended runtime: Python 3.6 with netCDF4, numba, numpy, e.g. Continuum.IO Anaconda
+"""
+import netCDF4 as nc4
+import numpy as np
+import numba
+
+def load(ncfilename, varname):
+    nc = nc4.Dataset(ncfilename)
+    nc.set_auto_scale(False)
+    ncv = nc[varname]
+    counts = ncv[:]
+    bit_depth = ncv.sensor_band_bit_depth
+    splat = counts.ravel()
+    return splat, bit_depth
+
+@numba.jit((numba.int16[:], numba.int32))
+def dist(splat, bit_depth):
+    z = np.zeros((2**bit_depth,), dtype=np.uint32)
+    for v in splat:
+        z[v] += 1
+    return z
+
+def hist(pngname, filename, distribution, extra=''):
+    import matplotlib
+    matplotlib.use('Agg')
+    import matplotlib.pyplot as plt
+    fig = plt.figure(figsize=(60,20))
+    ax = fig.add_subplot(111) #, axisbg='#EEEEEE')
+    # ax.grid(color='white', linestyle='solid')
+    ax.grid()
+    ax.plot(np.log(distribution+1), '.')
+    ax.set_xlabel('unscaled pixel value')
+    ax.set_ylabel('log(occurrences+1)')
+    ax.set_title('log(n+1) {} value distribution'.format(filename) + extra)
+    fig.savefig(pngname)
+
+if __name__=='__main__':
+    import sys, os
+    if len(sys.argv)==1:
+        print("parameters: netcdf-filename counts-variable-name optional-histogram optional-histogram-bins")
+    ncfilename, varname = sys.argv[1:3]
+    splat, bit_depth = load(ncfilename, varname)
+    distribution = dist(splat, bit_depth)
+    uniq = np.argwhere(distribution != 0).squeeze()
+    fn = os.path.split(ncfilename)[-1]
+    nvals = len(uniq)
+    print("{} of {} unique values in {}".format(nvals, 2**bit_depth, fn))
+    if len(sys.argv)>=4:
+        stem = sys.argv[3]
+        with open(stem + '.tsv', 'wt') as fp:
+            fp.write('\n'.join(['{}\t{}'.format(i,x) for i,x in enumerate(distribution)]) + '\n')
+        hist(stem + '.png', fn, distribution, ' {} unique values of {}'.format(nvals, 2**bit_depth))