ocssw V2020
proc_hico.py
Go to the documentation of this file.
1 #!/usr/bin/env python3
2 '''
3 Created on Nov 19, 2015
4 
5 @author: rhealy
6 '''
7 from netCDF4 import Dataset
8 
9 def proc_hico(ifile,earth_orient_file,leap_sec_file,boresight,nc_grp,calib_grp):
10  from hico.hicoLonLatHdr import hicoLonLatHdr
11  from hico.astreduc import frac2DegMinSec
12  import hico.bore_sight as bore
13  import hico.read_pos_vel_quat as rquat
14  import numpy as np
15  import os,sys
16  import struct
17  import hico.astreduc as astreduc
18 
19  lonHemi = {1:'E',0:'W'}
20  latHemi = {1:'N',0:'S'}
21  direction = {1:'ascending',0:'descending'}
22  byteorder = {'little':0, 'big':1}
23 
24  rad2Deg = 180.0/np.pi
25  deg2Rad = np.pi/180.0
26  theta_stowed=-60.
27 
28  try:
29  ocvarroot = os.environ['OCVARROOT']
30  except KeyError:
31  print("OCSSW environement variable OCVARROOT not set.")
32  print("Typically it is set to $OCSSW/run/var.")
33  sys.exit()
34 
35 
36  bs = np.zeros((3))
37  bs[0] = boresight[0]*deg2Rad
38  bs[1] = boresight[1]*deg2Rad
39  bs[2] = boresight[2]*deg2Rad
40  r_hico_bs = bore.bore_sight(bs)
41 
42  (quat_info, quat, pvq_data) = rquat.read_pos_vel_quat(ifile)
43  nsamples = np.int(quat_info['Requested number of pixels'])
44 # print('The CSV file is...' + quat_info['Source CSV file'])
45 # print('Theta = ' + quat_info['Theta (degrees from stowed position)'])
46 # print('The CSV file is...' + getattr(quat_info,'Source CSV file'))
47 # print('Distance Unit =' + getattr(quat, 'Distance Unit'))
48 
49  filein = "{}/hico/{}".format(ocvarroot,'nutation.dat')
50 # print(filein)
51  nut_data = astreduc.initReduc(filein)
52  astrodate = astreduc.HicoSec2Date(pvq_data['SecondsSinceEpoch'][0])
53 
54  yyyy = astrodate.yyyy % 100
55 
56  try:
57  pmx,pmy,dut1,loda = astreduc.getEDFData(earth_orient_file,yyyy,astrodate.mm,astrodate.dd)
58 # print(pmx,pmy,dut1,loda)
59  except ValueError:
60  print("Earth Data File does not have year={},month={},day={}".format(astrodate.yyyy,astrodate.mm,astrodate.dd))
61  sys.exit()
62 
63  try:
64 # print("Date={}/{}/{} {}:{}:{}".format(astrodate.yyyy,astrodate.mm,astrodate.dd,astrodate.hour,astrodate.min,astrodate.sec))
65  lsdat = astreduc.getLSData(leap_sec_file,astrodate.yyyy,astrodate.mm,astrodate.dd)
66 # print('lsdat=' + str(lsdat))
67  except ValueError:
68  print("Leap Sec File does not have year={},month={},day={}".format(astrodate.yyyy,astrodate.mm,astrodate.dd))
69  sys.exit()
70 
71 # Define the rotation from HICO frame of reference to ISS frame of reference
72 # NOTE: This assumes perfect alignment.
73 
74  #r_hico_to_iss = np.zeros((3,3))
75  #r_hico_to_iss[0][0] = -1.0
76  #r_hico_to_iss[1][1] = -1.0
77  #r_hico_to_iss[2][2] = 1.0
78  r_hico_to_iss = np.diagflat([-1]*2+[1])
79 # include bore sight offsets
80 
81  r_hico_to_iss = np.dot(r_hico_to_iss,r_hico_bs)
82 #!!!! BEGIN HICO pointing angles:
83 # I believe the 1/2 FOV is 3.461515 degrees (from Mike C?)
84 # A document by Bob Lucke states that the relationship is as follows;
85 # this was in MJM's HICO email folder in an email from Bill Snyder.
86 # UNITS ARE DEGREES
87 
88  if quat_info['ISS orientation'] == '+XVV':
89  imult = 1
90  else:
91  imult = -1
92  theta = np.zeros(nsamples)
93 
94  for i in range(0,nsamples):
95  frac = imult*(i-255.5)/255.5
96  theta[i] = ( -3.487 + 0.035 * frac**2)*frac
97 
98  theta = (theta + theta_stowed + float(quat_info['Theta (degrees from stowed position)']))*deg2Rad
99 # Now theta is in radians and accounts for angle from stowed position
100 # Now I need to convert to pointing angles in HICO's frame of reference.
101 # Note that HICO is in the YZ plane of HICO. The stowed angle (+ is right
102 # hand rule about X, measured from +Z) is in the +Y+Z quadrant, even
103 # for negative values of theta...
104 
105  svec = np.zeros((3,nsamples))
106 
107  svec[0,:] = 0.0
108  svec[1,:] = -np.sin(theta[:])
109  svec[2,:] = np.cos(theta[:])
110 
111  pm = astreduc.polarm(pmx,pmy)
112 
113  #print(pm)
114 
115  r_iss = np.array((3),dtype=float)
116  q_iss = np.array((4),dtype=float)
117  fout = quat_info['Expected Lon/Lat/View angle filename']
118  fheader = fout.split('.bil')[0] + '.hdr'
119  fbil=open(fout,"wb")
120  lines = pvq_data['SecondsSinceEpoch'].size
121 
122  #Setup the HICO header file using the hicoLonLatHdr class
123  hicoHeader = hicoLonLatHdr(filename=fheader, \
124  samples=nsamples, \
125  lines=lines, \
126  description='HICO geometry file, with calculated positions, and solar & view geometry, on the ground', \
127  sensor_type='HICO-ISS' \
128  )
129  hicoHeader.set_variable('data type',5)
130  hicoHeader.set_variable('interleave','bil')
131  hicoHeader.set_variable('wavelength units','{ degrees, degrees, degrees, degrees, degrees, degrees }')
132  hicoHeader.set_variable('band names','{ longitude, latitude, view zenith, view azimuth, sol zenith, sol azimuth }')
133  hicoHeader.set_variable('byte order',byteorder[sys.byteorder])
134  hicoHeader.set_variable('geometry_deltaUT1_s',dut1)
135  hicoHeader.set_variable('geometry_deltaT_s',lsdat)
136  hicoHeader.set_variable('geometry_xp_rad',pmx)
137  hicoHeader.set_variable('geometry_yp_rad',pmy)
138  hicoHeader.set_variable('geometry_LOD',loda)
139  hicoHeader.set_variable('geometry_bore_sight_params','{ ' + ', '.join( str(v) for v in bs) + ' }')
140 
141  for i in range(0,lines):
142 # for i in range(0,1):
143 
144  r_iss = np.multiply([pvq_data['ISSPOSX'][i],pvq_data['ISSPOSY'][i],pvq_data['ISSPOSZ'][i] ],0.3048e0) # convert to meeters
145  q_iss = [ pvq_data['ISSQS'][i],pvq_data['ISSQX'][i],pvq_data['ISSQY'][i],pvq_data['ISSQZ'][i] ]
146  rot_body = astreduc.quat_to_rot(q_iss)
147  astrodate = astreduc.HicoSec2Date(pvq_data['SecondsSinceEpoch'][i])
148  hico_jdate = astreduc.UT2time(astrodate.yyyy,astrodate.mm,astrodate.dd,astrodate.hour,astrodate.min,astrodate.sec,dut1,lsdat)
149  #print(hico_jdate.jdut1,hico_jdate.ttdb)
150  prec = astreduc.precession(hico_jdate.ttdb)
151  DeltaPsi, TrueEps, MeanEps, Omega, nut = astreduc.nutation(hico_jdate.ttdb,nut_data)
152  st,stdot,OmegaEarth = astreduc.sidereal(hico_jdate.jdut1,DeltaPsi,TrueEps,Omega,loda,2)
153 
154  t_eci_to_ecef=np.dot(pm,np.dot(st,np.dot(nut,prec)))
155  # for attitudes from HICO frame to ECEF
156  t_hico_to_ecef=np.dot(t_eci_to_ecef,np.dot(rot_body,r_hico_to_iss))
157  # The position does not depend on the attitude transformation
158  r_ecef=np.dot(t_eci_to_ecef,r_iss)
159  # Here, v_ecef are the 512 pointing vectors.
160  # v_ecef=matmul(t_hico_to_ecef,v_hico) ! v_hico=[3,NS]
161  v_ecef=np.dot(t_hico_to_ecef,svec)
162  llh,view_zen,view_az = astreduc.wgs84_intercept(r_ecef,v_ecef)
163 
164  sol_zen_az=rad2Deg*astreduc.solar_geometry(astrodate.yyyy,astrodate.mm,astrodate.dd,astrodate.hour,astrodate.min,astrodate.sec,llh[1,:],-llh[0,:])
165 
166  lon_out = np.array(llh[0,:])
167  lat_out = np.array(llh[1,:])
168  viewzen = np.array(view_zen,dtype='f8')
169  viewaz = np.array((360.0 +view_az) % 360.0,dtype='f8')
170  solarzen = np.array(sol_zen_az[0,:],dtype='f8')
171  solaraz = np.array(sol_zen_az[1,:],dtype='f8')
172  nc_grp['sensor_zenith'][i] = viewzen
173  nc_grp['sensor_azimuth'][i] = viewaz
174  nc_grp['solar_zenith'][i] = solarzen
175  nc_grp['solar_azimuth'][i] = solaraz
176  nc_grp['longitude'][i] = lon_out
177  nc_grp['latitude'][i] = lat_out
178  bilout = np.concatenate((lon_out,lat_out))
179  bilout = np.concatenate((bilout,viewzen))
180  bilout = np.concatenate((bilout,viewaz))
181  bilout = np.concatenate((bilout,solarzen))
182  bilout = np.concatenate((bilout,solaraz ))
183  myfmt='d'*len(bilout)
184  binout = struct.pack(myfmt,*bilout)
185  fbil.write(binout)
186 
187  # Below are some items for output header, grabbed near center of line
188  # Populate the HICO header for output to the header file
189  if i == lines/2-1:
190  hicoHeader.set_variable('date', '{ ' + str(astrodate.yyyy)+','+str(astrodate.mm)+','+str(astrodate.dd) + ' }')
191  hicoHeader.set_variable('time', '{ ' + str(astrodate.hour)+','+str(astrodate.min)+','+str(astrodate.sec) + ' }')
192  hicoHeader.set_variable('image_center_long','{ ' + ', '.join( str(v) for v in frac2DegMinSec(np.abs(lon_out[254]),0.,0.)) + ' }')
193  hicoHeader.set_variable('image_center_lat', '{ ' + ', '.join( str(v) for v in frac2DegMinSec(np.abs(lat_out[254]),0.,0.)) + ' }')
194  hicoHeader.set_variable('image_center_long_hem', lonHemi[(lon_out[254]>=0)])
195  hicoHeader.set_variable('image_center_lat_hem', latHemi[(lat_out[254]>=0)])
196  hicoHeader.set_variable('image_center_zenith_ang','{ ' + ', '.join( str(v) for v in frac2DegMinSec((viewzen[254]+viewzen[255])/2,0.,0.)) + ' }')
197  hicoHeader.set_variable('image_center_azimuth_ang','{ ' + ', '.join( str(v) for v in frac2DegMinSec((viewaz[254]+viewaz[255])/2,0.,0.)) + ' }')
198  hicoHeader.set_variable('geometry_ISS_Z_direction',direction[(pvq_data['ISSVELZ'][i]>=0)]) #! a bit crude, based on ECI velocity
199  r_iss_new = np.tile(r_iss[np.newaxis],1).T
200  llh = astreduc.ecef2latlon(r_iss_new)
201  hicoHeader.set_variable('sensor_altitude',llh[2,0]/1000.)
202 
203  fbil.close()
204 
205  # Create the history for output into the header file using the information from the
206  # pos_vel_quat input csv file and then write out the header file.
207 
208  hicoHeader.createHistory(quat_info)
209  hicoHeader.writeHicoENVIHeaderFile()
210 
211  calib_grp.hico_orientation_from_quaternion = quat_info['ISS orientation']
212 
213 if __name__ == '__main__':
214  import argparse
215 
216  parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter, description='''\
217  Generate a HICO geometry file.
218 
219  ''', add_help=True)
220  parser.add_argument('-ifile', nargs=1, type=str, help=' iss*pos_vel_quat.csv input file (Must be a CSV file) ')
221  parser.add_argument('-lfile', nargs=1, type=str, help=' leapsec file ')
222  parser.add_argument('-efile', nargs=1, type=str, help=' earth-orient file ')
223  parser.add_argument('-boresight', nargs=3, type=float, default=([0, 0, 0]), help=('Process Bore Sight Parameters'))
224 
225  args = parser.parse_args('-ifile iss.2013067.0308.063527.L0.12933.20130308205743.hico_pos_vel_quat.csv -efile /home/rhealy/ocsswn/run/var/hico/finals.data -lfile /home/rhealy/ocsswn/run/var/modis/leapsec.dat -boresight -0.9957 0.0268 -0.0128'.split())
226 
227  ifile = ''.join(args.ifile)
228  earth_orient_file = ''.join(args.efile)
229  leap_sec_file = ''.join(args.lfile)
230  boresight = args.boresight
231  rootgrp = Dataset('test_hico.nc','w')
232  rootgrp.createDimension('scan_lines', 2000)
233  rootgrp.createDimension('samples', 512)
234  #
235  # L2gen is expecting these variables under group "navigation" for l2gen
236  #
237  navigation = rootgrp.createGroup("navigation")
238 
239  senz = navigation.createVariable('sensor_zenith','f4',('scan_lines','samples',))
240  solz = navigation.createVariable('solar_zenith','f4',('scan_lines','samples',))
241  sena = navigation.createVariable('sensor_azimuth','f4',('scan_lines','samples',))
242  sola = navigation.createVariable('solar_azimuth','f4',('scan_lines','samples',))
243  lon = navigation.createVariable('longitude','f4',('scan_lines','samples',))
244  lat = navigation.createVariable('latitude','f4',('scan_lines','samples',))
245  senz.units = 'degrees'
246  senz.valid_min = 0
247  senz.valid_max = 180
248  senz.long_name = 'sensor zenith'
249  solz.units = 'degrees'
250  solz.valid_min = 0
251  solz.long_name = 'solar zenith'
252  solz.valid_max = 180
253  sena.units = 'degrees'
254  sena.valid_min = -180
255  sena.valid_max = 180
256  sena.long_name = 'sensor azimuth'
257  sola.units = 'degrees'
258  sola.valid_min = -180
259  sola.valid_max = 180
260  sola.long_name = 'solar azimuth'
261  lon.units = 'degrees'
262  lon.valid_min = -180
263  lon.valid_max = 180
264  lon.long_name = 'longitude'
265  lat.units = 'degrees'
266  lat.valid_min = -180
267  lat.valid_max = 180
268  lat.long_name = 'latitude'
269  # orientation needs to be set in group = "metadata/HICO/Calibration for l2gen"
270  metadata = rootgrp.createGroup("metadata")
271  hico = metadata.createGroup("HICO")
272  calibration = hico.createGroup("Calibration")
273  proc_hico(ifile,earth_orient_file,leap_sec_file,boresight,navigation,calibration)
274 
275  rootgrp.close()
def frac2DegMinSec(fracd, fracm, fracs)
Definition: astreduc.py:1135
def proc_hico(ifile, earth_orient_file, leap_sec_file, boresight, nc_grp, calib_grp)
Definition: proc_hico.py:9