Ich habe ein Problem beim Plotten eines RGB-Bildes mit Pythons Grundkartenmodul mit Breiten- und Längengraddaten. Jetzt kann ich die Plots machen, die ich will, aber das Problem ist, wie langsam es ist, da es in der Lage ist, einzelne Kanaldaten viel schneller als die RGB-Daten zu plotten, und im Allgemeinen RGB-Bilder selbst zu zeichnen schnell. Da ich lat/lon-Daten habe, werden die Dinge kompliziert. Ich habe die Lösung für dieses Problem ausgecheckt:So zeichnen Sie Geolocated-RGB-Daten schneller mit Python-Grundkarte
How to plot an irregular spaced RGB image using python and basemap?
das ist, wie ich da, wo ich jetzt bin. Es geht im Wesentlichen um das folgende Problem. Wenn Sie in der Grundkarte die Methode pcolormesh verwenden, um RGB-Daten zu plotten, müssen Sie einen colorTuple-Parameter definieren, der die RGB-Daten Punkt für Punkt abbildet. Da die Array-Größe in der Größenordnung von 2000x1000 liegt, dauert dies eine Weile. Ein Ausschnitt von dem, was ich spreche, ist unten zu sehen (volle Arbeits Code weiter unten):
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)
Wenn nur ein Kanal Plotten, kann er die Karte in etwa 10 Sekunden machen oder so. Beim Plotten der RGB-Daten kann es 3-4 Minuten dauern. Angesichts der Tatsache, dass es nur 3-mal so viele Daten gibt, glaube ich, dass es einen besseren Weg geben muss, zumal das Zeichnen von RGB-Daten bei rechteckigen Bildern genauso schnell gehen kann wie die Daten eines einzelnen Kanals.
Meine Fragen sind also: Gibt es eine Möglichkeit, diese Berechnung schneller zu machen, entweder mit anderen Plot-Modulen (Bokeh zum Beispiel) oder durch Ändern der Farbabbildung in irgendeiner Weise? Ich habe versucht, imshow mit sorgfältig ausgewählten Kartengrenzen zu verwenden, aber da es nur das Bild im vollen Umfang der Karte streckt, ist dies nicht wirklich gut genug für eine genaue Zuordnung der Daten.
Unten ist eine abgespeckte Version von meinem Code, der für ein Beispiel mit der richtigen Module arbeiten:
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)