Cómo trazar datos RGB geolocalizados más rápido con el mapa base de Python

Tengo un problema al trazar una imagen RGB utilizando el módulo de mapa base de Python con datos de latitud y longitud. Ahora, puedo hacer los gráficos que quiero, pero el problema es lo lento que es, ya que puede trazar datos de un solo canal mucho más rápido que los datos RGB, y en general, trazar imágenes RGB por su cuenta también es rápido. Como tengo datos de lat / lon, ahí es donde las cosas se complican. He comprobado la solución a este problema:

¿Cómo trazar una imagen RGB con espacios irregulares usando python y el mapa base?

que es como llegué a donde estoy ahora. Básicamente se trata de la siguiente cuestión. Al usar el método pcolormesh en el mapa base, para trazar datos RGB, debe definir un parámetro colorTuple que mapeará los datos RGB punto por punto. Dado que el tamaño de la matriz es del orden de 2000×1000, esto toma un tiempo para hacerlo. A continuación se muestra un fragmento de lo que estoy hablando (código de trabajo completo más abajo):

 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) 

Al trazar un solo canal, puede hacer el mapa en aproximadamente 10 segundos aproximadamente. Al trazar los datos RGB, puede tomar de 3 a 4 minutos. Dado que solo hay 3 veces más datos, creo que debe haber una mejor manera, especialmente porque la representación de datos RGB puede ir tan rápido como los datos de un canal cuando se crean imágenes rectangulares.

Entonces, mi pregunta es: ¿Hay alguna forma de hacer este cálculo más rápido, ya sea con otros módulos de trazado (por ejemplo, Bokeh) o cambiando el mapeo de colores de alguna manera? He intentado usar imshow con los límites del mapa cuidadosamente seleccionados, pero como solo se extiende la imagen en toda la extensión del mapa, esto no es lo suficientemente bueno para un mapeo preciso de los datos.

A continuación se muestra una versión simplificada de mi código que funcionará para un ejemplo con los módulos correctos:

 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) 

Resolví este problema agregando el 1.0 al valor de cada parte del colorTuple para convertirlo en una matriz RGBA. pcolormesh función pcolormesh y descubrí que estaba llamando al convertidor de color para convertir el RGB a una matriz RGBA 4 veces diferentes, lo que me tomaba unos 50 segundos cada vez. Si le da una matriz RGBA para que comience, omitirá esto y producirá la ttwig en un período de tiempo razonable. La línea de código adicional que se agregó se ve a continuación:

 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)