3
我遇到一個問題，使用具有緯度和經度數據的Python底圖模塊繪製RGB圖像。現在，我可以製作出我想要的圖，但問題在於它的速度有多慢，因爲它能夠繪製比RGB數據快得多的單通道數據，而且一般情況下，單獨繪製RGB圖像也是快速。由於我有緯度/經度數據，這就是事情變得複雜的地方。我已經簽出瞭解決這個問題：如何使用Python底圖更快地繪製地理定位的RGB數據
How to plot an irregular spaced RGB image using python and basemap?
這是我得到了我現在在哪裏。它基本上歸結爲以下問題。在底圖中使用pcolormesh方法時，要繪製RGB數據，您必須定義colorTuple參數，該參數將逐點映射RGB數據。由於數組大小約爲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數據的速度可能與單通道數據一樣快。
所以，我的問題是：有沒有辦法讓其他繪圖模塊（例如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)