From ada75c9e98357d93e73cee19bbe9fd66d0c8e4ac Mon Sep 17 00:00:00 2001 From: RKGarcia <rayg@ssec.wisc.edu> Date: Thu, 6 Apr 2017 02:03:18 +0000 Subject: [PATCH] rad_count_dist with plotting utility --- goesr/rad_count_dist.py | 58 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 goesr/rad_count_dist.py diff --git a/goesr/rad_count_dist.py b/goesr/rad_count_dist.py new file mode 100644 index 0000000..04570a5 --- /dev/null +++ b/goesr/rad_count_dist.py @@ -0,0 +1,58 @@ +#!/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)) -- GitLab