python画rgb图_如何使用Python底图更快地绘制地理定位的RGB数据

frompyhdf.SDimportSD,SDCimportnumpyasnpfromscipy.interpolateimportinterp1dimportmatplotlib.pyplotaspltfrommpl_toolkits.basemapimportBasemapdefget_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()returnattr_outdefget_hdf_dataset(infile,dataset):f=SD(infile,SDC.READ)data=f.select(dataset)[:]f.end()returndataclassmake_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.0rgb[rgb<0]=0.0lin=np.array([0,30,60,120,190,255])/255.0nonlin=np.array([0,110,160,210,240,255])/255.0scale=interp1d(lin,nonlin,kind='quadratic')self.img=scale(rgb)defplot_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()defplot_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)ifone_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 Fastmake_rgb(data_file).plot_image()# Also Fast, takes about 10 secondsmake_rgb(data_file).plot_geo(geo_file,one_channel=True)# Much slower, takes several minutesmake_rgb(data_file).plot_geo(geo_file)