2016-12-30 4 views
3

緯度と経度のデータを持つPythonのBasemapモジュールを使用してRGBイメージをプロットする際に問題があります。さて、私は必要なプロットを作ることができましたが、問題はどれほど遅いのですか?単一チャンネルのデータをRGBデータよりもはるかに高速にプロットすることができ、一般的にRGB画像をプロットすることも可能です速い。私には緯度/経度データがあるので、状況は複雑になります。私は、この問題への解決策をチェックアウトしました:Pythonのベースマップを使用してジオロケーション付きRGBデータをより速くプロットする方法

How to plot an irregular spaced RGB image using python and basemap?

私は今午前どこに得た方法です。それは基本的に以下の問題になります。 pcolormeshメソッドをベースマップで使用する場合、RGBデータをプロットするには、RGBデータをポイントごとにマッピングするcolorTupleパラメータを定義する必要があります。配列のサイズは2000x1000程度なので、これにはしばらく時間がかかります。以下に見られる私が話しているのスニペット(完全な作業コード、さらにダウン):ただ一つのチャネルをプロットすると

if one_channel: 
    m.pcolormesh(lons, lats, img[:,:,0], latlon=True) 
else: 
    # This is the part that is slow, but I don't know how to 
    # accurately plot the data otherwise. 

    mesh_rgb = img[:, :-1, :] 
    colorTuple = mesh_rgb.reshape((mesh_rgb.shape[0] * mesh_rgb.shape[1]), 3) 

    # What you put in for the image doesn't matter because of the color mapping 
    m.pcolormesh(lons, lats, img[:,:,0], latlon=True,color=colorTuple) 

、それは約10秒ほどでマップを作ることができます。 RGBデータをプロットすると、3〜4分かかります。 3倍のデータしかないことを考えれば、特に、RGBデータをプロットすると、矩形の画像を作成するときに1つのチャンネルデータと同じ速さで処理できるため、より良い方法が必要であると感じています。

私の質問は次のとおりです。他のプロットモジュール(Bokehなど)を使って、あるいは色のマッピングを変更することで、この計算を高速化する方法はありますか?私はimshowを慎重に選択したマップ境界で試してみましたが、画像をマップの全範囲に伸ばしているだけなので、データの正確なマッピングには十分ではありません。以下は

正しいモジュールと例のために動作します私のコードのストリップダウンバージョンです:

from pyhdf.SD import SD,SDC 
import numpy as np 
from scipy.interpolate import interp1d 
import matplotlib.pyplot as plt 
from mpl_toolkits.basemap import Basemap 

def get_hdf_attr(infile,dataset,attr): 

    f = SD(infile,SDC.READ) 
    data = f.select(dataset) 
    index = data.attr(attr).index() 
    attr_out = data.attr(index).get() 
    f.end() 

    return attr_out 

def get_hdf_dataset(infile,dataset): 

    f = SD(infile,SDC.READ) 
    data = f.select(dataset)[:] 
    f.end() 

    return data 

class make_rgb: 

    def __init__(self,file_name): 

     sds_250 = get_hdf_dataset(file_name, 'EV_250_Aggr1km_RefSB') 
     scales_250 = get_hdf_attr(file_name, 'EV_250_Aggr1km_RefSB', 'reflectance_scales') 
     offsets_250 = get_hdf_attr(file_name, 'EV_250_Aggr1km_RefSB', 'reflectance_offsets') 

     sds_500 = get_hdf_dataset(file_name, 'EV_500_Aggr1km_RefSB') 
     scales_500 = get_hdf_attr(file_name, 'EV_500_Aggr1km_RefSB', 'reflectance_scales') 
     offsets_500 = get_hdf_attr(file_name, 'EV_500_Aggr1km_RefSB', 'reflectance_offsets') 

     data_shape = sds_250.shape 

     along_track = data_shape[1] 
     cross_track = data_shape[2] 

     rgb = np.zeros((along_track, cross_track, 3)) 

     rgb[:, :, 0] = (sds_250[0, :, :] - offsets_250[0]) * scales_250[0] 
     rgb[:, :, 1] = (sds_500[1, :, :] - offsets_500[1]) * scales_500[1] 
     rgb[:, :, 2] = (sds_500[0, :, :] - offsets_500[0]) * scales_500[0] 

     rgb[rgb > 1] = 1.0 
     rgb[rgb < 0] = 0.0 

     lin = np.array([0, 30, 60, 120, 190, 255])/255.0 
     nonlin = np.array([0, 110, 160, 210, 240, 255])/255.0 

     scale = interp1d(lin, nonlin, kind='quadratic') 

     self.img = scale(rgb) 

    def plot_image(self): 

     fig = plt.figure(figsize=(10, 10)) 
     ax = fig.add_subplot(111) 

     ax.set_yticks([]) 
     ax.set_xticks([]) 
     plt.imshow(self.img, interpolation='nearest') 
     plt.show() 

    def plot_geo(self,geo_file,one_channel=False): 

     fig = plt.figure(figsize=(10, 10)) 
     ax = fig.add_subplot(111) 

     lats = get_hdf_dataset(geo_file, 0) 
     lons = get_hdf_dataset(geo_file, 1) 

     lat_0 = np.mean(lats) 
     lat_range = [np.min(lats), np.max(lats)] 

     lon_0 = np.mean(lons) 
     lon_range = [np.min(lons), np.max(lons)] 

     map_kwargs = dict(projection='cass', resolution='l', 
          llcrnrlat=lat_range[0], urcrnrlat=lat_range[1], 
          llcrnrlon=lon_range[0], urcrnrlon=lon_range[1], 
          lat_0=lat_0, lon_0=lon_0) 

     m = Basemap(**map_kwargs) 

     if one_channel: 
      m.pcolormesh(lons, lats, self.img[:,:,0], latlon=True) 
     else: 
      # This is the part that is slow, but I don't know how to 
      # accurately plot the data otherwise. 
      mesh_rgb = self.img[:, :-1, :] 
      colorTuple = mesh_rgb.reshape((mesh_rgb.shape[0] * mesh_rgb.shape[1]), 3) 
      m.pcolormesh(lons, lats, self.img[:,:,0], latlon=True,color=colorTuple) 

     m.drawcoastlines() 
     m.drawcountries() 

     plt.show() 

if __name__ == '__main__': 

    # https://ladsweb.nascom.nasa.gov/archive/allData/6/MOD021KM/2015/183/ 
    data_file = 'MOD021KM.A2015183.1005.006.2015183195350.hdf' 

    # https://ladsweb.nascom.nasa.gov/archive/allData/6/MOD03/2015/183/ 
    geo_file = 'MOD03.A2015183.1005.006.2015183192656.hdf' 

    # Very Fast 
    make_rgb(data_file).plot_image() 

    # Also Fast, takes about 10 seconds 
    make_rgb(data_file).plot_geo(geo_file,one_channel=True) 

    # Much slower, takes several minutes 
    make_rgb(data_file).plot_geo(geo_file) 

答えて

3

私はオンにcolorTupleのすべての部分の値に1.0を追加することで、この問題を解決しましたそれをRGBA配列に変換します。私はpcolormesh関数を調べ、カラーコンバーターを呼び出してRGBAをRGBA配列に4回ずつ変換し、毎回約50秒かかると判明しました。開始するRGBA配列を与えると、これをバイパスして妥当な時間枠でプロットを生成します。追加された追加のコード行は次のとおりです。

if one_channel: 
    m.pcolormesh(lons, lats, img[:,:,0], latlon=True) 
else: 
    mesh_rgb = img[:, :-1, :] 
    colorTuple = mesh_rgb.reshape((mesh_rgb.shape[0] * mesh_rgb.shape[1]), 3) 

    # ADDED THIS LINE 
    colorTuple = np.insert(colorTuple,3,1.0,axis=1) 

    # What you put in for the image doesn't matter because of the color mapping 
    m.pcolormesh(lons, lats, img[:,:,0], latlon=True,color=colorTuple) 
関連する問題