How to merge multiple MODIS swaths into one plot in python?

Question

I want to mosaic/merge multiple swaths of the MODIS dataset (MOD06_L2) using python. I used the example (http://hdfeos.org/zoo/MORE/LAADS/MOD/MOD04_L2_merge.py) to read multiple files and merge. But I am getting an error while doing so, how to correct this error?

I would like to know is there any better way than this, to merge/mosaic MODIS HDF files into one?

import os
import glob                                                                 
import matplotlib as mpl
import matplotlib.pyplot as plt
# import cartopy.crs as ccrs
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np

# The first file in 3 swath files.
FILE_NAME = 'MOD06_L2.A2017126.0655.061.2019226193408.hdf'
GEO_FILE_NAME ='MOD06_L2.A2017126.0655.061.2019226193408.hdf'
DATAFIELD_NAME = 'Brightness_Temperature'

from pyhdf.SD import SD, SDC

i = 0

for file in list(glob.glob('MOD06*.hdf')):
    reader = open(file)
    hdf = SD(file, SDC.READ)
    # Read dataset.
    data2D = hdf.select(DATAFIELD_NAME)
    data = data2D[:,:].astype(np.double)
    hdf_geo = SD(GEO_FILE_NAME, SDC.READ)
    # Read geolocation dataset.
    lat = hdf_geo.select('Latitude')
    latitude = lat[:,:]
    lon = hdf_geo.select('Longitude')
    longitude = lon[:,:]
    # Retrieve attributes.
    attrs = data2D.attributes(full=1)
    lna=attrs["long_name"]
    long_name = lna[0]
    aoa=attrs["add_offset"]
    add_offset = aoa[0]
    fva=attrs["_FillValue"]
    _FillValue = fva[0]
    sfa=attrs["scale_factor"]
    scale_factor = sfa[0]
    vra=attrs["valid_range"]
    valid_min = vra[0][0]
    valid_max = vra[0][1]        
    ua=attrs["units"]
    units = ua[0]
    invalid = np.logical_or(data > valid_max,data < valid_min)
    invalid = np.logical_or(invalid, data == _FillValue)
    data[invalid] = np.nan
    data = (data - add_offset) * scale_factor 
    datam = np.ma.masked_array(data, np.isnan(data))
    if i == 0 :
        data_m = datam
        latitude_m = latitude
        longitude_m = longitude
    else:
        data_m = np.vstack([data_m, datam])
        latitude_m = np.vstack([latitude_m, latitude])
        longitude_m = np.vstack([longitude_m, longitude])
    i = i + 1
m = Basemap(projection='cyl', resolution='l',
                llcrnrlat=-90, urcrnrlat=90,
                llcrnrlon=-180, urcrnrlon=180)
m.drawcoastlines(linewidth=0.5)
m.drawparallels(np.arange(-90, 91, 45))
m.drawmeridians(np.arange(-180, 180, 45), labels=[True,False,False,True])
sc = m.scatter(longitude_m, latitude_m, c=data_m, s=0.1, cmap=plt.cm.jet,
               edgecolors=None, linewidth=0)
cb = m.colorbar()
cb.set_label(units)


# Put title using the first file.
basename = os.path.basename(FILE_NAME)
plt.title('{0}\n{1}'.format(basename, DATAFIELD_NAME))
fig = plt.gcf()

# Save image.
pngfile = "{0}.py.png".format(basename)
fig.savefig(pngfile)

It showing an error

ValueError: 'c' argument has 4604040 elements, which is inconsistent with 'x' and 'y' with size 657720.