Q:

GDAL couldn't create raster with equal shape

I am trying to carry out pixel by pixel correlation for each image in a time series vs another 1d array of data. To do that, the arrays coming out of the raster files must be identical in shape. However, when I use the code below, I get different shapes. 
def change_raster_extent(extent, in_file, arrays):
    out_tif = in_file.replace('.tif', '_ext.tif')
    ds = gdal.Open(in_file)
    gdal.Translate(out_tif, ds, outputBounds=extent)

    band = ds.GetRasterBand(1)
    # dis_arr = np.reshape(discharge, (21, 0)).T
    x_array = np.array(band.ReadAsArray())
    print x_array.shape
    arrays.append(x_array)
    return out_tif

extent = [-98.05392182  17.40680802 -80.43304074  31.46484375]
arrays = []
for in_file in in_files:
    change_raster_extent(extent, in_file, arrays)

It outputs:
(971, 1217)
(982, 1231)
(985, 1234)
(985, 1235)
(991, 1242)
(995, 1247)
(993, 1245)
(967, 1213)
(992, 1244)
(978, 1226)
(986, 1236)
(972, 1219)

Edit: I realized the possible cause of this issue is because these images are MODIS aqua satellite data and their pixel positions don't match. So even though their geographic extent is same, the pixels don match. I suspect, re-sampling the images could fix the issue. I will try that and see if it fixes it. 
How could I make the shape to be identical for the output array?

A:

I had to use this solution that mainly uses gdalbuildvrt. The function below aligns rasters based on the first raster in a list. 
def align_rasters(rasters, output_suffix):
    """
    Aligns list of rasters to have the same resolution and 
    cell size for pixel based calculations. 
    :param rasters: List of raster paths.
    :type rasters: List
    :param output_suffix: The output aligned rasters files suffix with extension.
    :type output_suffix: String
    :return: True if the process runs and False if the data couldn't be read. 
    :rtype: Boolean
    """
    command = ["gdalbuildvrt", "-te"]
    hDataset = gdal.Open(rasters[0], gdal.GA_ReadOnly)
    if hDataset is None:
        return False
    adfGeoTransform = hDataset.GetGeoTransform(can_return_null=True)
    for tif_file in rasters:

        vrt_file = tif_file.replace('.tif', '.vrt')

        if adfGeoTransform is not None:
            dfGeoXUL = adfGeoTransform[0]
            dfGeoYUL = adfGeoTransform[3]
            dfGeoXLR = adfGeoTransform[0] + adfGeoTransform[1] * hDataset.RasterXSize + \
                       adfGeoTransform[2] * hDataset.RasterYSize
            dfGeoYLR = adfGeoTransform[3] + adfGeoTransform[4] * hDataset.RasterXSize + \
                       adfGeoTransform[5] * hDataset.RasterYSize
            xres = str(abs(adfGeoTransform[1]))
            yres = str(abs(adfGeoTransform[5]))
            subprocess.call(command + [str(dfGeoXUL), str(dfGeoYLR), str(dfGeoXLR),
                                       str(dfGeoYUL), "-q", "-tr", xres, yres,
                                       vrt_file, tif_file])
            output_file = tif_file.replace('.tif', output_suffix)
            cmd = 'gdal_translate -q {} {}'.format(vrt_file, output_file)
            subprocess.call(cmd, shell=True)
            os.remove(vrt_file)
    return True

The above function requires list of rasters to be aligned and a suffix for the output.
Example usage: align_rasters(['path/to/raster1.tif', 'path/to/raster1.tif'], '_aligned.tif')

