Rudimentary satellite predictions and line identification

imgstat.py

@@ -8,6 +8,7 @@ from astropy.coordinates import SkyCoord, AltAz, EarthLocation
 from astropy.time import Time
 import configparser
 import argparse
+import os
 
 # Read commandline options
 conf_parser = argparse.ArgumentParser(description='Plot image statistics')
@@ -19,6 +20,10 @@ conf_parser.add_argument("-i", "--input",
                          help="Specify file to be processed. If no file" +
                          " is specified ./imgstat.csv will be used.",
                          metavar='FILE', default="./imgstat.csv")
+conf_parser.add_argument("-d", "--directory",
+                         help="Specify directory of observations. If no" +
+                         " directory is specified parent will be used.",
+                         metavar='DIR', dest='file_dir', default=".")
 conf_parser.add_argument("-o", "--output",
                          help="Specify output file. Default is 'imgstat.png'",
                          metavar='FILE', default="./imgstat.png")
@@ -32,6 +37,9 @@ if args.conf_file:
 else:
     cfg.read('configuration.ini')
 
+# Move to processing directory
+os.chdir(args.file_dir)
+    
 table = ascii.read(args.input, format="csv")
 
 t = Time(table['mjd'], format="mjd", scale="utc")

process.py

@@ -5,6 +5,8 @@ import numpy as np
 from stvid.stio import fourframe
 from stvid.stars import generate_star_catalog
 from stvid.astrometry import calibrate_from_reference
+from stvid.satellite import generate_satellite_predictions
+from stvid.satellite import find_hough3d_lines
 import astropy.units as u
 from astropy.utils.exceptions import AstropyWarning
 from astropy.coordinates import EarthLocation
@@ -45,11 +47,14 @@ if __name__ == "__main__":
                         lon=cfg.getfloat('Common', 'observer_lon')*u.deg,
                         height=cfg.getfloat('Common', 'observer_el')*u.m)
 
+    # Move to processing directory
+    os.chdir(args.file_dir)
+    
     # Get files
-    files = sorted(glob.glob(os.path.join(args.file_dir, "2*.fits")))
+    files = sorted(glob.glob("2*.fits"))
 
     # Statistics file
-    fstat = open(os.path.join(args.file_dir, "imgstat.csv"), "w")
+    fstat = open("imgstat.csv", "w")
     fstat.write("fname,mjd,ra,de,rmsx,rmsy,mean,std,nstars,nused\n")
 
     # Loop over files
@@ -58,10 +63,15 @@ if __name__ == "__main__":
         pix_catalog = generate_star_catalog(fname)
 
         # Calibrate astrometry
-        calibrate_from_reference(fname,
-                                 os.path.join(args.file_dir, "test.fits"),
+        calibrate_from_reference(fname, "test.fits",
                                  pix_catalog)
 
+        # Generate satellite predictions
+        generate_satellite_predictions(fname)
+
+        # Extract lines with 3D Hough transform
+        ids=find_hough3d_lines(fname)
+        
         # Stars available and used
         nused = np.sum(pix_catalog.flag == 1)
         nstars = pix_catalog.nstars

stvid/satellite.py

@@ -0,0 +1,65 @@
+#!/usr/bin/env python
+from __future__ import print_function
+import os
+import subprocess
+from stvid.stio import fourframe
+from stvid.stio import satid
+from stvid.stio import observation
+from tempfile import mkstemp
+
+def generate_satellite_predictions(fname):
+    # Format command
+    command = "satid %s %s.png/png"%(fname, fname)
+
+    # Run command
+    output = subprocess.check_output(command, shell=True,
+                                     stderr=subprocess.STDOUT)
+
+    return
+
+def find_hough3d_lines(fname,ntrkmin=20,dr=8):
+    # Read four frame
+    ff=fourframe(fname)
+
+    # Mask frame
+    ff.mask(10,10,5,5)
+
+    # Compute selection mask
+    x,y,z,t,sig=ff.selection_mask(5.0,40.0)
+    
+    # Save points to temporary file
+    f,fpath=mkstemp()
+    with open(fpath,"w") as f:
+        for i in range(len(t)):
+            f.write("%f,%f,%f\n"%(x[i],y[i],z[i]))
+    f.close()
+
+    # Run 3D Hough line-finding algorithm
+    command="hough3dlines -dx %d -minvotes %d %s"%(dr,ntrkmin,fpath)
+    output=subprocess.check_output(command,shell=True,stderr=subprocess.STDOUT)
+
+    # Remove file
+    os.remove(fpath)
+
+    # Clean output (a bit cluncky)
+    cleaned_output = output.replace("npoints=", "")
+    cleaned_output = cleaned_output.replace(", a=(", " ")
+    cleaned_output = cleaned_output.replace("), b=("," ")
+    cleaned_output = cleaned_output.replace(")", "")
+    cleaned_output = cleaned_output.replace(",", " ")
+    print(cleaned_output)
+    # Generate identifications
+    lines=[]
+    s=cleaned_output.split()
+    for i in range(len(s)//7):
+        x0,y0,z0=float(s[1+7*i]),float(s[2+7*i]),float(s[3+7*i])
+        dx,dy,dz=float(s[4+7*i]),float(s[5+7*i]),float(s[6+7*i])
+        xmin=x0-z0*dx/dz
+        xmax=x0+(ff.nz-z0)*dx/dz
+        ymin=y0-z0*dy/dz
+        ymax=y0+(ff.nz-z0)*dy/dz
+        line="%.23s %8.3f %8.3f %8.3f %8.3f %8.5f  %s unidentified sunlit"%(ff.nfd,xmin,ymin,xmax,ymax,ff.texp,99999)
+        lines.append(line)
+
+    return [satid(line) for line in lines]
+

stvid/stars.py

@@ -4,7 +4,6 @@ import os
 import subprocess
 import numpy as np
 
-
 class pixel_catalog:
     """Pixel catalog"""
 
@@ -32,7 +31,7 @@ class pixel_catalog:
 
 def generate_star_catalog(fname):
     # Skip if file already exists
-    if not os.path.exists(fname+".cat"):
+    if not os.path.exists(os.path.join(fname, ".cat")):
         # Get sextractor location
         env = os.getenv("ST_DATADIR")
 

stvid/stio.py

@@ -4,7 +4,7 @@ import numpy as np
 from astropy.io import fits
 from astropy.time import Time
 from astropy import wcs
-
+from scipy import ndimage
 
 class observation:
     """Satellite observation"""
@@ -93,6 +93,9 @@ class fourframe:
             # Read image planes
             self.zavg, self.zstd, self.zmax, self.znum = hdu[0].data
 
+            # Generate sigma frame
+            self.zsig=(self.zmax-self.zavg)/(self.zstd+1e-9)
+            
             # Frame properties
             self.ny, self.nx = self.zavg.shape
             self.nz = hdu[0].header['NFRAMES']
@@ -141,6 +144,46 @@ class fourframe:
         self.w.wcs.ctype = self.ctype
         self.w.wcs.set_pv([(2, 1, 45.0)])
 
+
+    def mask(self,xmin,xmax,ymin,ymax):
+        x,y=np.meshgrid(np.arange(self.nx),np.arange(self.ny))
+        c=(x>=xmin) & (x<=self.nx-xmax) & (y>=ymin) & (y<=self.ny-ymax)
+        self.mask=np.ones_like(self.zavg)
+        self.mask[~c]=0.0
+        self.zavg*=self.mask
+        self.zstd*=self.mask
+        self.zmax*=self.mask
+        self.znum*=self.mask
+        self.zsig*=self.mask
+
+    def selection_mask(self,sigma,zstd):
+        """Create a selection mask"""
+        c1=ndimage.uniform_filter(self.znum,3,mode='constant')
+        c2=ndimage.uniform_filter(self.znum*self.znum,3,mode='constant')
+        z=np.sqrt(c2-c1*c1)
+
+        # Standard deviation mask
+        c=z<zstd
+        m1=np.zeros_like(self.zavg)
+        m1[c]=1.0
+
+        # Sigma mask
+        c=self.zsig<sigma
+        m2=np.zeros_like(self.zavg)
+        m2[~c]=1.0
+        self.zsel=m1*m2
+
+        # Generate points
+        c=self.zsel==1.0
+        xm,ym=np.meshgrid(np.arange(self.nx),np.arange(self.ny))
+        x,y=np.ravel(xm[c]),np.ravel(ym[c])
+        inum=np.ravel(self.znum[c]).astype('int')
+        sig=np.ravel(self.zsig[c])
+        t=np.array([self.dt[i] for i in inum])
+
+        return x,y,inum,t,sig
+
+        
     def significant_pixels_along_track(self, sigma, x0, y0,
                                        dxdt, dydt, rmin=10.0):
         """Extract significant pixels along a track"""