Python OCO_Matrix.data[0,band_idx+1]方法代码示例

# 需要导入模块: from OCO_Matrix import OCO_Matrix [as 别名]
# 或者: from OCO_Matrix.OCO_Matrix import data[0,band_idx+1] [as 别名]
def create_scene_dispersion_file(sounding_id, latitude, sza_r, saz_r, time_struct, aband_data, dispersion_coefs, apriori_out_file, l1b_build_id=None, index_scheme=None):
   """
   Purpose: Given a single frame of gosat data, estimate the spectral shift in band 1P.
            This is accomplished by using the strong solar line at 12985.163 wavenumbers.
            This is in the a-band continuum and is really the only local feature there.
   
            Note: This routine fits for a global shift, which INCLUDES the instrument
            doppler shift.  If this is not desired, the user must subtract off
            the instrument doppler shift.
   
   Original Author: Chris Odell
   Converted from IDL
   """

   if type(sounding_id) is str:
      pol_name = sounding_id[-1]
      if pol_name not in ACOS_File.GOSAT_POL_ORDER:
         # If not P or S then use averaged offset
         pol_name = 'averaged'
   else:
      raise Exception('Can not determine polarization name from non string ;Asounding id: %s' % sounding_id)

   wn_s0 = ABAND_SOLAR_LINE
   wn_s0 += ABAND_ILS_OFFSET[pol_name] # account for different ILS offset P vs S

   # Calculate the Speed of Earth center towards sun
   frac = (time_struct.tm_hour + time_struct.tm_min/60. + time_struct.tm_sec/3600.)/24. # fraction of a day
   doy = time_struct.tm_yday + frac + 0.5 # why 0.5? well cause L2 code does that
   V_cen = 497.2 * math.sin((doy-4.1)/365.25*2*math.pi) # simple approximate model, good to ~ 10 m/s.

   # Calculate the rotational component of the solar doppler velocity
   geo_lat = math.atan(math.tan(math.radians(latitude))/(1.+CON))
   Rloc = 1000.0 * RADIUS/math.sqrt(1.+ CON*math.sin(geo_lat)**2)
   V_rot = -EARTH_ROT_SPD * Rloc * math.sin(sza_r) * math.cos(saz_r-math.radians(90.)) * math.cos(geo_lat)

   solar_doppler_velocity = V_cen + V_rot
  
   # Modify position of strong solar line to include solar doppler shift
   wn_s = wn_s0 * (1.0e0 - solar_doppler_velocity/constants.c) # takes into account the solar doppler shift

   # Not used but could be good for debugging?
   #solar_shift = wn_s - wn_s0

   # Create aband dispersion array, 1 indexed
   abo2_dcoef_0, abo2_dcoef_1 = dispersion_coefs[0, :]
   aband_disp = numpy.arange(1, len(aband_data)+1)*abo2_dcoef_1 + abo2_dcoef_0

   # (2) CONSTRUCT THE X, Y FUNCTION TO FIT
   #
   # Use intersection of the points above and below range threshold
   # Make sure to sort the resulting indexes since set.intersection
   # does not guarantee anything about ordering
   w_1 = numpy.where(aband_disp >= ABAND_DISP_SEARCH_RANGE[0])
   w_2 = numpy.where(aband_disp <= ABAND_DISP_SEARCH_RANGE[1])
   w = sorted(list(set.intersection(set(w_1[0]), set(w_2[0]))))

   if len(w) == 0:
      raise ValueError('Could not find any points in the range %s in the aband dispersion with range %s for sounding %s' % (ABAND_DISP_SEARCH_RANGE, (min(aband_disp), max(aband_disp)), sounding_id))

   x = aband_disp[w] - wn_s

   mxmeas = max(aband_data[w]) # maximum value
   m2 = max(mxmeas - aband_data[w])
   y = (mxmeas-aband_data[w])/m2 # this should look like a gaussian

   w = numpy.where(y < 0.1)
   y = y - numpy.mean(y[w]) # subtract off the "continuum"
   pos = numpy.argmax(y) # pos = index of the maximum value of y
   m = len(y)
   fg = -x[pos]  # first-guess value of spectral shift

   # (3) PERFORM A SIMPLE FIT ASSUMING A PERFECTLY LINEAR MODEL
   #     MODEL IS A GAUSSIAN WITH SIGMA=0.2 CM^-1
   #     FIT PARAMETERS ARE AMPLITUDE AND CENTER OF THE GAUSSIAN
   K = numpy.zeros((m, 2), dtype=float) # will hold jacobian
   sig2 = 0.2e0**2
   ygauss = numpy.exp(-(x+fg)**2/ (2.0*sig2))
   K[:, 0] = ygauss
   K[:, 1] = -ygauss / sig2 * (x+fg)
   Kt = numpy.transpose(K)

   # Note that KtK is a 2x2 symmetric matrix and can be inverted analytically if desired.
   # form matrix (Kt K)^{-1} Kt
   delta = numpy.dot(numpy.linalg.inv(numpy.dot(Kt, K)), (numpy.dot(Kt, (y-ygauss))))
   fit = [1.0, fg] + delta

   wn_shift = fit[1]

   # Debugging to make sure fitting works
   if (0):
      from matplotlib import pyplot
      print sounding_id, aband_disp[0], wn_shift, aband_disp[0]-aband_disp[1]

      pyplot.cla()
      pyplot.plot(x, y) # measured data
      yfit = fit[0]* numpy.exp(-(x+fit[1])**2/ (2.0*sig2))
      pyplot.plot(x,yfit) # fit modeled data

      pyplot.legend(('Measured Solar Line', 'Fitted Solar Line'), loc='lower left')
      
#.........这里部分代码省略.........