simpler, robust API

ConvertDMC2h5.py

@@ -22,7 +22,8 @@ def dmclooper(p):
     cmosinit = {'firstrawind':p.cmos[0],'lastrawind':p.cmos[1]}
 
     params = {'kineticsec':p.kineticsec,'rotccw':p.rotccw,'transpose':p.transpose,
-              'flipud':p.flipud,'fliplr':p.fliplr,'fire':p.fire,'sensorloc':p.loc}
+              'flipud':p.flipud,'fliplr':p.fliplr,'fire':p.fire,'sensorloc':p.loc,
+              'cmdlog':' '.join(argv)}
 
     infn = Path(p.infile).expanduser()
     if infn.is_file():
@@ -45,7 +46,7 @@ def dmclooper(p):
         rawImgData,rawind,finf = goRead(f, p.pix,p.bin,p.frames,p.ut1,
                                     p.kineticsec,p.startutc,cmosinit,p.verbose,ofn,p.headerbytes)
 #%% convert
-        vid2h5(None, finf['ut1'], rawind, None, ofn, params, argv)
+        vid2h5(None, finf['ut1'], rawind, None, ofn, params)
 #%% optional plot
         if p.movie:
             plots(rawImgData,rawind,finf)

histutils/__init__.py

@@ -6,6 +6,7 @@
 import numpy as np
 import logging
 from struct import pack,unpack
+from typing import Tuple
 # NOTE: need to use int64 since Numpy thru 1.11 defaults to int32 on Windows for dtype=int, and we need int64 for large files
 
 
@@ -25,7 +26,7 @@ def write_quota(outbytes:int, outfn:Path) -> int:
         return freeout
 
 
-def sixteen2eight(I,Clim):
+def sixteen2eight(I:np.ndarray, Clim:tuple) -> np.ndarray:
     """
     scipy.misc.bytescale had bugs
 
@@ -39,7 +40,8 @@ def sixteen2eight(I,Clim):
     Q *= 255 # stretch to [0,255] as a float
     return Q.round().astype(np.uint8) # convert to uint8
 
-def normframe(I,Clim):
+
+def normframe(I:np.ndarray, Clim:tuple) -> np.ndarray:
     """
     inputs:
     -------
@@ -51,8 +53,7 @@ def normframe(I,Clim):
     return (I.astype(np.float32).clip(Vmin, Vmax) - Vmin) / (Vmax - Vmin) #stretch to [0,1]
 
 
-
-def getRawInd(fn:Path, BytesPerImage:int, nHeadBytes:int, Nmetadata:int):
+def getRawInd(fn:Path, BytesPerImage:int, nHeadBytes:int, Nmetadata:int) -> Tuple[int,int]:
     assert isinstance(Nmetadata,int)
     if Nmetadata<1: #no header, only raw images
         fileSizeBytes = fn.stat().st_size
@@ -107,6 +108,7 @@ def req2frame(req, N:int=0):
 
     return frame
 
+
 def dir2fn(ofn,ifn,suffix) -> Path:
     """
     ofn = filename or output directory, to create filename based on ifn
@@ -133,6 +135,7 @@ def dir2fn(ofn,ifn,suffix) -> Path:
 
     return ofn
 
+
 def splitconf(conf,key,i=None,dtype=float,fallback=None,sep=','):
     if conf is None:
         return fallback
@@ -177,18 +180,21 @@ def splitconf(conf,key,i=None,dtype=float,fallback=None,sep=','):
         else:
             return fallback
 # %% HDF5
-
 def setupimgh5(f, Nframetotal:int, Nrow:int, Ncol:int, dtype=np.uint16,
-               writemode='r+', key='/rawimg',cmdlog:str=''):
+               writemode='r+', key='/rawimg',cmdlog:str=None):
     """
     f: HDF5 handle (or filename)
 
     h: HDF5 dataset handle
     """
     if isinstance(f, (str,Path)):  # assume new HDF5 file wanted
+        if f.is_dir():
+            raise IOError('must provide specific HDF5 filename, not just a directory')
+
         print('creating', f)
-        with h5py.File(f, writemode, libver='latest') as F:
-            return setupimgh5(F,Nframetotal,Nrow,Ncol,dtype,writemode,key)
+        with h5py.File(f, writemode) as F:
+            setupimgh5(F,Nframetotal,Nrow,Ncol,dtype,writemode,key)
+
     elif isinstance(f, h5py.File):
         h = f.create_dataset(key,
                  shape =  (Nframetotal,Nrow,Ncol),
@@ -210,9 +216,8 @@ def setupimgh5(f, Nframetotal:int, Nrow:int, Ncol:int, dtype=np.uint16,
     else:
         raise TypeError(f'{type(f)} is not correct, must be filename or h5py.File HDF5 file handle')
 
-    return h
 
-def vid2h5(data, ut1, rawind, ticks, outfn, P, cmdlog='', i:int=0, Nfile:int=1, det=None,tstart=None):
+def vid2h5(data:np.ndarray, ut1, rawind, ticks, outfn:Path, P:dict, i:int=0, Nfile:int=1, det=None,tstart=None):
     assert outfn,'must provide a filename to write'
 
     outfn = Path(outfn).expanduser()
@@ -253,7 +258,7 @@ def vid2h5(data, ut1, rawind, ticks, outfn, P, cmdlog='', i:int=0, Nfile:int=1,
                     break
             ind = slice(None)
 
-        with h5py.File(outfn, writemode, libver='latest') as f:
+        with h5py.File(outfn, writemode) as f:
             if data is not None:
                 if 'rawimg' not in f:  # first run
                     setupimgh5(f, N, data.shape[1], data.shape[2])
@@ -355,8 +360,9 @@ def vid2h5(data, ut1, rawind, ticks, outfn, P, cmdlog='', i:int=0, Nfile:int=1,
         savemat(outfn,matdata,oned_as='column')
     else:
         raise ValueError(f'what kind of file is {outfn}')
+
 
-def imgwriteincr(fn, imgs, imgslice):
+def imgwriteincr(fn:Path, imgs, imgslice):
     """
     writes HDF5 huge image files in increments
     """
@@ -367,9 +373,9 @@ def imgwriteincr(fn, imgs, imgslice):
     if isinstance(fn, Path):
         assert fn.suffix == '.h5','Expecting to write .h5 file' # avoid accidental overwriting of source file due to misspecified command line
 
-        with h5py.File(fn, 'r+', libver='latest') as f:
+        with h5py.File(fn, 'r+') as f:
             f['/rawimg'][imgslice,:,:] = imgs
     elif isinstance(fn, h5py.File):
         f['/rawimg'][imgslice,:,:] = imgs
     else:
-        raise TypeError(f'fn must be Path or h5py.File instead of {type(fn)}')
+        raise TypeError(f'{fn} must be Path or h5py.File instead of {type(fn)}')

histutils/camclass.py

@@ -6,19 +6,23 @@
 from scipy.signal import savgol_filter
 from numpy.random import poisson
 import h5py
+from typing import List
 #
 from . import splitconf
 from pymap3d import azel2radec, aer2ecef
-from pymap3d.haversine import angledist,angledist_meeus
-from dascutils.readDASCfits import readDASC
+from pymap3d.haversine import anglesep
+import dascutils.io as dio
 from themisasi.fov import mergefov
 from .plots import plotnear_rc, plotlsq_rc
 
 class Cam: #use this like an advanced version of Matlab struct
     """
     simple mode via try except attributeerror
     """
-    def __init__(self,sim,cp,name,zmax=None,xreq=None,makeplot=[],calfn=None,verbose=0):
+    def __init__(self,sim, cp:dict, name:str,
+                 zmax=None,xreq=None,makeplot:List[str]=[],
+                 calfn:Path=None,verbose:int=0):
+
         self.verbose = verbose
 
         try:
@@ -49,9 +53,9 @@ def __init__(self,sim,cp,name,zmax=None,xreq=None,makeplot=[],calfn=None,verbose
             if splitconf(cp,'plotMaxVal',ci) is not None:
                 self.clim[1] =  splitconf(cp,'plotMaxVal',ci)
 
-            _,(self.az,self.el),self.lla,_ = readDASC(self.fn[0],
-                                                cp['azcalfn'].split(',')[ci],
-                                                cp['elcalfn'].split(',')[ci])
+            _,(self.az,self.el),self.lla,_ = dio.load(self.fn[0],
+                                                        cp['azcalfn'].split(',')[ci],
+                                                        cp['elcalfn'].split(',')[ci])
 
             if 'realvid' in makeplot:
                 if 'h5' in makeplot:
@@ -71,9 +75,15 @@ def __init__(self,sim,cp,name,zmax=None,xreq=None,makeplot=[],calfn=None,verbose
             return
 
         elif self.usecam:
-            self.name = int(name)
+            try: # integer name
+                self.name = int(name)
+            except ValueError: # non-integer name
+                self.name = name
 #%%
-        self.ncutpix = int(cp['nCutPix'].split(',')[ci])
+        if isinstance(cp['nCutPix'],str):
+            self.ncutpix = int(cp['nCutPix'].split(',')[ci])
+        else: # int
+            self.ncutpix = cp['nCutPix']
 
         self.Bincl =  splitconf(cp,'Bincl',ci)
         self.Bdecl =  splitconf(cp,'Bdecl',ci)
@@ -144,11 +154,16 @@ def __init__(self,sim,cp,name,zmax=None,xreq=None,makeplot=[],calfn=None,verbose
 
         # data file name
         if sim.realdata and self.usecam:
-            #.strip() needed in case of multi-line ini
-            self.fn = (sim.realdatapath / cp['fn'].split(',')[ci].strip()).expanduser()
+            if isinstance(cp['fn'],str):
+                #.strip() needed in case of multi-line ini
+                self.fn = (sim.realdatapath / cp['fn'].split(',')[ci].strip()).expanduser()
+            elif isinstance(cp['fn'],Path):
+                self.fn = sim.realdatapath / cp['fn']
+            else:
+                raise ValueError('your camera filename {cp["fn"]} not found')
 
             if not self.fn.suffix == '.h5':
-                raise TypeError('I can only work with HDF5 files, not FITS. Use ConvertSolisFITSh5 if you have FITS')
+                raise OSError('I can only work with HDF5 files, not FITS. Use demutils/ConvertSolisFITSh5 if you have FITS')
 
             assert self.fn.is_file(),'{} does not exist'.format(self.fn)
 
@@ -414,21 +429,19 @@ def findLSQ(self,nearrow,nearcol,odir):
         assert len(radecMagzen) == 2
         logging.info('mag. zen. ra,dec {}'.format(radecMagzen))
 
-        if False: #astropy
-            angledist_deg = angledist(      radecMagzen[0],radecMagzen[1], self.ra[cutrow, cutcol], self.dec[cutrow, cutcol])
-        else: #meeus
-            angledist_deg = angledist_meeus(radecMagzen[0],radecMagzen[1], self.ra[cutrow, cutcol], self.dec[cutrow, cutcol])
+        anglesep_deg = anglesep(radecMagzen[0],radecMagzen[1], self.ra[cutrow, cutcol], self.dec[cutrow, cutcol])
+
 #%% assemble angular distances
-        self.angle_deg,self.angleMagzenind = self.sky2beam(angledist_deg)
+        self.angle_deg,self.angleMagzenind = self.sky2beam(anglesep_deg)
 
         self.cutrow = cutrow
         self.cutcol = cutcol
 #%% meager self-check of result
         if self.arbfov is not None:
             expect_diffang = self.arbfov / self.ncutpix
             diffang = np.diff(self.angle_deg)
-            diffoutlier = max(abs(expect_diffang-diffang.min()),
-                              abs(expect_diffang-diffang.max()))
+           # diffoutlier = max(abs(expect_diffang-diffang.min()),
+           #                   abs(expect_diffang-diffang.max()))
             assert_allclose(expect_diffang, diffang.mean(), rtol=0.01),'large bias in camera angle vector detected'
 #            assert diffoutlier < expect_diffang,'large jump in camera angle vector detected' #TODO arbitrary
 

histutils/findnearest.py

@@ -12,6 +12,10 @@ def findClosestAzel(az,el,azpts,elpts):
     """
     assert az.ndim     == 2
     assert az.shape    == el.shape
+
+    azpts = np.atleast_1d(azpts)
+    elpts = np.atleast_1d(elpts)
+
     assert azpts.shape == elpts.shape
 
     az = np.ma.masked_invalid(az)

histutils/get1Dcut.py

@@ -1,10 +1,13 @@
+from pathlib import Path
 from numpy import logspace
 import h5py
+from typing import List
 #
 from pymap3d import ecef2aer
 from .plots import plotLOSecef
+from .camclass import Cam
 
-def get1Dcut(cam,odir,verbose):
+def get1Dcut(cam:List[Cam], odir:Path=None, verbose:bool=False) -> List[Cam]:
     """
     i.   get az/el of each pixel (rotated/transposed as appropriate)
     ii.  get cartesian ECEF of each pixel end, a point outside the grid (to create rays to check intersections with grid)
@@ -37,7 +40,7 @@ def get1Dcut(cam,odir,verbose):
     if verbose and odir:
         dbgfn = odir / 'debugLSQ.h5'
         print('writing', dbgfn)
-        with h5py.File(str(dbgfn),'w',libver='latest') as f:
+        with h5py.File(dbgfn,'w') as f:
             for C in cam:
                 f[f'/cam{C.name}/cutrow'] = C.cutrow
                 f[f'/cam{C.name}/cutcol'] = C.cutcol