0

I need to create an RGB figure based on an HDF5 file, which consists of latitude layer, longitude layer, and different band data layers. However, the ax.contourf() function only receive 2-dimensions(one band) data instead of 3-dimensions data (RGB,3 bands); ax.pcolormesh and ax.imshow() can not receive latitude/longitude parameters. Now, I can achieve my goal by using ax.inshow() function with the GeoTIFF file converted from the HDF5 file. But, how can I achieve that by using an HDF5 file instead of Geotiff?

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from matplotlib.colors import LogNorm
import matplotlib.colors as colors
import pandas as pd
import glob,matplotlib
import numpy as np
import georaster
from PIL import Image,ImageOps

class stretch: 

    def hist_calc(img, ratio):
        img[img < 0]=np.nan
        
        hist, bins = np.histogram(img, range(int(np.nanmin(img)),int(np.nanmax(img))))
        total_pixels =np.sum(hist)
        
        min_index = int(ratio * total_pixels)
        max_index = int((1 - ratio) * total_pixels)
        min_gray = 0
        max_gray = 0
        
        sum1 = 0
        for i in range(hist.__len__()):
            sum1 = sum1 + hist[i]
            if sum1 > min_index:
                min_gray = i
                break
        
        sum2 = 0
        for i in range(hist.__len__()):
            sum2= sum2 + hist[i]
            if sum2 > max_index:
                max_gray = i
                break
        return min_gray, max_gray

    def linearStretch(array, ratio):
        
        old_min, old_max = stretch.hist_calc(array, ratio)
    
        array[array < old_min]=old_min
        array[array > old_max]=old_max
#         print('old min = %d,old max = %d' % (np.nanmin(array), np.nanmax(array)))        
        array_min, array_max = np.nanmin(array), np.nanmax(array)
        return ((array - array_min)/(array_max - array_min))


if __name__=='__main__':
    files=glob.glob(r'H:\calibration\DL\HY1114\H1B_OPER_OCT_*20080305T_024104696_geo_ROI_RC_cloudmask.img')
    for i,img in enumerate(files):
        red =georaster.SingleBandRaster(img,band=8)
        green =georaster.SingleBandRaster(img,band=7)
        blue =georaster.SingleBandRaster(img,band=5)
        
        fig=plt.figure(11,figsize=(2.9,3.1))
        ax=plt.axes([0.05,0.05,0.94,0.94])
        b_map = Basemap(resolution='l', area_thresh=100, projection='cyl',\
                        llcrnrlon = red.extent[0]-3, llcrnrlat = red.extent[2]-3, urcrnrlon = red.extent[1]+3, urcrnrlat = red.extent[3]+3)  
        matplotlib.rcParams["font.family"] = "Times New Roman"  
        csfont = {'fontname':'Times New Roman','fontsize':10}
        parallels = [red.extent[2],red.extent[3]]
        b_map.drawparallels(parallels,labels=[True,False,False,False],rotation='vertical',**csfont)
        meridians = [red.extent[0],red.extent[1]]
       
        b_map.drawmeridians(meridians,labels=[False,False,False,True],yoffset=-2,**csfont)
        b_map.drawcoastlines(linewidth=0.2,color='red')

        nl,ns=(red.r).shape[0],(red.r).shape[1]
        raster=np.full([nl,ns,3], np.nan)
        raster[:, :,0]=stretch.linearStretch(red.r, 0.05)
        raster[:, :,1]=stretch.linearStretch(green.r,  0.05)
        raster[:, :,2]=stretch.linearStretch(blue.r,  0.05)

        rgb = np.dstack((raster[:, :,0], raster[:, :,1], raster[:, :,2]))
        plt.imshow(rgb, extent=red.extent)#,
        figname=os.path.dirname(img)+'/'+ os.path.basename(img)[0:-3]+'jpg' 
        plt.savefig(figname,dpi=300)
        plt.show()

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.