2015-11-03 57 views
1

我試圖繪製數據以檢查我的代碼,我將所得繪圖與已經用Matlab生成的繪圖進行比較。然而,我遇到了幾個問題:在Python中繪製RINEX

  1. 一般來說,RINEX文件的解析工作,並且數據表示的一般模式看起來類似於Matlab腳本繪製的。但是,在數據放大時,數據會有小的偏差,例如在使用較小的時間序列時,例如繪製特殊的2小時而非24小時。在Matlab中,可以看到這個小的差異,並應用多項式擬合。但是對於Python圖(下面顯示的第一張圖),這兩小時期間的曲線看起來「平滑」,並且完全沒有偏離,就像在Matlab腳本中看到的那樣(第二個圖顯示藍色線作爲數據,相對於polyfit的紅線,因此,藍線在x = 9.4時表現出輕微的差異)。假定Matlab腳本是正確的,因爲這種偏差是由於地震活動暫時破壞了電離層。請參閱下面的情節:

Python plot of data

Matlab plot of data

第三個情節是在Matlab,其中這簡直是polyfit減去實時數據。

因此,不清楚該數據是如何繪製在Python腳本的座標軸上的,因爲數據看起來很平滑?也不如果我的代碼是錯誤的(見下文)並以某種方式「平滑」出的數據以某種方式:

#Calculating by looping through 
      for sv in range(32): 
        sat = self.obs_data_chunks_dataframe[sv, :] 
        #print "sat.index_{0}: {1}".format(sv+1, sat.index) 
        phi1 = sat['L1'] * LAMBDA_1 #Change units of L1 to meters 
        phi2 = sat['L2'] * LAMBDA_2 #Change units of L2 to meters 
        pr1 = sat['P1'] 
        pr2 = sat['P2'] 


        #CALCULATION: teqc Calculation 
        iono_teqc = COEFF * (pr2 - pr1)/1000000 #divide to make values smaller (tbc) 
        print "iono_teqc_{0}: {1}".format(sv+1, iono_teqc) 


        #PLOTTING 
        #Plotting of the data 
        plt.plot(sat.index, iono_teqc, label=‘teqc’) 
        plt.xlabel('Time (UTC)') 
        plt.ylabel('Ionosphere Delay (meters)') 
        plt.title("Ionosphere Delay on {0} for Satellite {1}.".format(self.date, sv+1)) 
        plt.legend() 
        ax = plt.gca() 
        ax.ticklabel_format(useOffset=False) 
        plt.grid() 

        if sys.platform.startswith('win'): 
         plt.savefig(winpath + '\Figure_SV{0}'.format(sv+1)) 
        elif sys.platform.startswith('darwin'): 
         plt.savefig(macpath + 'Figure_SV{0}'.format(sv+1)) 

        plt.close() 
  • 從點1繼,多項式擬合代碼下面不按照我想要的方式運行,所以我在這裏忽略了一些東西。我認爲這與在x,y軸上使用的數據有關,但不能精確地確定是什麼。有人會知道我哪裏錯了嗎?

     #Zoomed in plots 
         if sv == 19: 
          #Plotting of the data 
          plt.plot(sat.index, iono_teqc, label=‘teqc’) #sat.index to plot for time in UTC 
          plt.xlim(8, 10) 
          plt.xlabel('Time (UTC)') 
          plt.ylabel('Ionosphere Delay (meters)') 
          plt.title("Ionosphere Delay on {0} for Satellite {1}.".format(self.date, sv+1)) 
          plt.legend() 
          ax = plt.gca() 
          ax.ticklabel_format(useOffset=False) 
          plt.grid() 
    
          #Polynomial fitting 
          coefficients = np.polyfit(sat.index, iono_teqc, 2) 
          plt.plot(coefficients) 
    
          if sys.platform.startswith('win'): 
           #os.path.join(winpath, 'Figure_SV{0}'.format(sv+1)) 
           plt.savefig(winpath + '\Zoom_SV{0}'.format(sv+1)) 
          elif sys.platform.startswith('darwin'): 
           plt.savefig(macpath + 'Zoom_SV{0}'.format(sv+1)) 
    
          plt.close() 
    
  • 我的RINEX文件包含32顆衛星。然而試圖生成的曲線對於所有32時,我得到:

    IndexError:指數31超出範圍爲軸線0與尺寸31

  • 變更下方至31個解決了這個部分中,只有代碼不包括第32顆衛星。我想也爲衛星32繪圖。解析和格式化數據的功能如下:

    def read_obs(self, RINEXfile, n_sat, sat_map): 
        obs = np.empty((TOTAL_SATS, len(self.obs_types)), dtype=np.float64) * np.NaN 
        lli = np.zeros((TOTAL_SATS, len(self.obs_types)), dtype=np.uint8) 
        signal_strength = np.zeros((TOTAL_SATS, len(self.obs_types)), dtype=np.uint8) 
    
        for i in range(n_sat): 
         # Join together observations for a single satellite if split across lines. 
         obs_line = ''.join(padline(RINEXfile.readline()[:-1], 16) for _ in range((len(self.obs_types) + 4)/5)) 
         #obs_line = ''.join(padline(RINEXfile.readline()[:-1], 16) for _ in range(2)) 
         #while obs_line 
    
    
    
         for j in range(len(self.obs_types)): 
          obs_record = obs_line[16*j:16*(j+1)] 
          obs[sat_map[i], j] = floatornan(obs_record[0:14]) 
          lli[sat_map[i], j] = digitorzero(obs_record[14:15]) 
          signal_strength[sat_map[i], j] = digitorzero(obs_record[15:16]) 
    
        return obs, lli, signal_strength 
    
    
    def read_data_chunk(self, RINEXfile, CHUNK_SIZE = 10000): 
        obss = np.empty((CHUNK_SIZE, TOTAL_SATS, len(self.obs_types)), dtype=np.float64) * np.NaN 
        llis = np.zeros((CHUNK_SIZE, TOTAL_SATS, len(self.obs_types)), dtype=np.uint8) 
        signal_strengths = np.zeros((CHUNK_SIZE, TOTAL_SATS, len(self.obs_types)), dtype=np.uint8) 
        epochs = np.zeros(CHUNK_SIZE, dtype='datetime64[us]') 
        flags = np.zeros(CHUNK_SIZE, dtype=np.uint8) 
    
        i = 0 #ggfrfg 
        while True: 
         hdr = self.read_epoch_header(RINEXfile) 
         if hdr is None: 
          break 
         epoch_time, flags[i], sats = hdr 
         #epochs[i] = np.datetime64(epoch_time) 
         epochs[i] = epoch_time 
         sat_map = np.ones(len(sats)) * -1 
         for n, sat in enumerate(sats): 
          if sat[0] == 'G': 
           sat_map[n] = int(sat[1:]) - 1 
         obss[i], llis[i], signal_strengths[i] = self.read_obs(RINEXfile, len(sats), sat_map) 
         i += 1 
         if i >= CHUNK_SIZE: 
          break 
    
        return obss[:i], llis[:i], signal_strengths[:i], epochs[:i], flags[:i] 
    
    
    def read_data(self, RINEXfile): 
        obs_data_chunks = [] 
    
        while True: 
         obss, _, _, epochs, _ = self.read_data_chunk(RINEXfile) 
         epochs = epochs.astype(np.int64) 
         epochs = np.divide(epochs, float(3600.000)) 
    
         if obss.shape[0] == 0: 
          break 
    
         obs_data_chunks.append(pd.Panel(
          np.rollaxis(obss, 1, 0), 
          items=['G%02d' % d for d in range(1, 33)], 
          major_axis=epochs, 
          minor_axis=self.obs_types 
         ).dropna(axis=0, how='all').dropna(axis=2, how='all')) 
    
         self.obs_data_chunks_dataframe = obs_data_chunks[0] 
    

    有什麼建議嗎?

    乾杯,pymat。

    回答

    0

    我設法解決了Qu1,因爲這是一個轉換問題,我的計算被忽略了,另外兩個點卻是開放的...