Python OCO_Matrix.header['solar_doppler_velocity']方法代码示例

# 需要导入模块: from OCO_Matrix import OCO_Matrix [as 别名]
# 或者: from OCO_Matrix.OCO_Matrix import header['solar_doppler_velocity'] [as 别名]

#.........这里部分代码省略.........
   # 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')
      
      pyplot.savefig(os.path.join(os.path.basename(apriori_out_file), 'disp_fit_%s.png' % sounding_id), format='png')


   disp_obj = OCO_Matrix()

   disp_obj.data = numpy.zeros((dispersion_coefs.shape[1], dispersion_coefs.shape[0]+1), dtype=float)

   # First column is index of parameter
   disp_obj.data[:, 0] = numpy.arange(1,dispersion_coefs.shape[1]+1)

   # Create dispersion file scaling non aband values according to this method:
   # Where below wn_spc is the spacing, ie second coefficient
   # wco2 offset : delta_2 = 0.48 * ( delta_1 + wn_spc) - wn_spc)
   # sco2_offset : delta_3 = 0.38 * ( delta_1 + wn_spc) - wn_spc)

   disp_offsets = []
   
   for band_idx, band_scaling in enumerate(DISPERSION_OFFSET_SCALING):
      # Assuming 2 coefficients for describing dispersion for GOSAT
      # Since we need both to do modifcation to offset
      band_wn_start, band_spacing = dispersion_coefs[band_idx, :]

      # Use a different equation for l1b builds before v02.04.06
      if (l1b_build_id != None and (0,0,0) > l1b_build_id < (2,4,6)) or (index_scheme != None and index_scheme == 0):
         # For 0 based indexing in L1B file
         disp_obj.header['scaling_equation'] = '"0 based indexing"'
         curr_offset = DISPERSION_OFFSET_SCALING[band_idx] * (wn_shift - band_spacing) + band_spacing
      else:
         # For 1 based indexing in L1B file
         disp_obj.header['scaling_equation'] = '"1 based indexing"'
         curr_offset = DISPERSION_OFFSET_SCALING[band_idx] * wn_shift
      
      disp_offsets.append(curr_offset)

      # First column is index of parameter
      disp_obj.data[0, band_idx+1] = band_wn_start + curr_offset
      disp_obj.data[1, band_idx+1] = band_spacing

   disp_obj.header['offsets'] = disp_offsets
   disp_obj.header['solar_doppler_velocity'] = solar_doppler_velocity
   disp_obj.labels = DISPERSION_LABELS
   disp_obj.file_id = 'Dispersion apriori for: %s' % sounding_id
   disp_obj.write(apriori_out_file)