2016-03-10 27 views
1

我對python非常陌生,目前我正試圖將CALIPSO反向散射數據繪製爲高度和位置的函數,而不是高度和時間。我不確定從哪裏開始,因爲我對此很陌生。 下面是我與繪製緯度/經度的CALIPSO反向散射數據函數

import numpy as np 
    import matplotlib as mpl 
    import matplotlib.pyplot as plt 
    from ccplot.hdf import HDF 
    from ccplot.algorithms import interp2d_12 
    import ccplot.utils 

    filename = '/CAL_LID_L1-ValStage1-V3-02.2012-12-01T04-24-22ZN.hdf' 
    name = 'Total_Attenuated_Backscatter_532' 
    label = 'Total Attenuated Backscatter 532nm (km$^{-1}$ sr$^{-1}$)' 
    colormap = '/ccplot/cmap/calipso-backscatter.cmap' 
    x1 = 0 
    x2 = 1000 
    h1 = 0 # km 
    h2 = 20 # km 
    nz = 500 # Number of pixels in the vertical. 

    if __name__ == '__main__': 
with HDF(filename) as product: 
    # Import datasets. 
    time = product['Profile_UTC_Time'][x1:x2, 0] 
    height = product['metadata']['Lidar_Data_Altitudes'] 
    dataset = product[name][x1:x2] 

    # Convert time to datetime. 
    time = np.array([ccplot.utils.calipso_time2dt(t) for t in time]) 

    # Mask missing values. 
    dataset = np.ma.masked_equal(dataset, -9999) 

    # Interpolate data on a regular grid. 
    X = np.arange(x1, x2, dtype=np.float32) 
    Z, null = np.meshgrid(height, X) 
    data = interp2d_12(
     dataset[::], 
     X.astype(np.float32), 
     Z.astype(np.float32), 
     x1, x2, x2 - x1, 
     h2, h1, nz, 
    ) 

    # Import colormap. 
    cmap = ccplot.utils.cmap(colormap) 
    cm = mpl.colors.ListedColormap(cmap['colors']/255.0) 
    cm.set_under(cmap['under']/255.0) 
    cm.set_over(cmap['over']/255.0) 
    cm.set_bad(cmap['bad']/255.0) 
    norm = mpl.colors.BoundaryNorm(cmap['bounds'], cm.N) 

    # Plot figure. 
    fig = plt.figure(figsize=(12, 6)) 
    TIME_FORMAT = '%e %b %Y %H:%M:%S UTC' 
    im = plt.imshow(
     data.T, 
     extent=(mpl.dates.date2num(time[0]), mpl.dates.date2num(time[-1]), h1, h2), 
     cmap=cm, 
     norm=norm, 
     aspect='auto', 
     interpolation='nearest', 
    ) 
    ax = im.axes 
    ax.set_title('CALIPSO %s - %s' % (
     time[0].strftime(TIME_FORMAT), 
     time[-1].strftime(TIME_FORMAT) 
    )) 
    ax.set_xlabel('Time') 
    ax.set_ylabel('Altitude (km)') 
    ax.xaxis.set_major_locator(mpl.dates.AutoDateLocator()) 
    ax.xaxis.set_major_formatter(mpl.dates.DateFormatter('%H:%M:%S')) 
    cbar = plt.colorbar(
     extend='both', 
     use_gridspec=True 
    ) 
    cbar.set_label(label) 
    fig.tight_layout() 
    plt.savefig('calipso-plot.png') 
    plt.show() 

回答

-1
  1. 工作使用獲得的緯度/經度信息的代碼:

    緯度=產物[ '緯度'] [X1:X2,0]
    經度=產品[ '經度'] [X1:X2,0]

  2. 下#劇情圖:

    IM = plt.imshow(..., 程度=(latitude.min(),latitude.max(),H1,H2), ...)

  3. 註釋出:

    ax.xaxis.set_major_locator(...) 斧.xaxis.set_major_formatter(...)

,並替換:

ax.xaxis.set_majot_locator(mpl.ticker.LinearLocator()) 

玩得開心!