Align two rasters so that data overlaps by geographical location
    Usage:
    (alignedraster_o, alignedraster_a, geot_a) = AlignRasters(raster, alignraster, how=np.mean)

    where:
        raster: string with location of raster to be aligned

        alignraster: string with location of raster to which raster will be aligned

        how: function used to aggregate cells (if the rasters have different sizes)

    It is assumed that both rasters have the same size

def align_rasters(raster, alignraster, how=np.ma.mean, cxsize=None, cysize=None, masked=False):
    '''
    Align two rasters so that data overlaps by geographical location
    Usage:
    (alignedraster_o, alignedraster_a, geot_a) = AlignRasters(raster, alignraster, how=np.mean)

    where:
        raster: string with location of raster to be aligned

        alignraster: string with location of raster to which raster will be aligned

        how: function used to aggregate cells (if the rasters have different sizes)

    It is assumed that both rasters have the same size
    '''
    ndv1, xsize1, ysize1, geot1, projection1, datatype1 = get_geo_info(raster)
    ndv2, xsize2, ysize2, geot2, projection2, datatype2 = get_geo_info(alignraster)
    if projection1.ExportToMICoordSys() == projection2.ExportToMICoordSys():
        blocksize = (np.round(geot2[1]/geot1[1]).astype(int), np.round(geot2[-1]/geot1[-1]).astype(int))
        mraster = gdalnumeric.LoadFile(raster)
        mraster = np.ma.masked_array(mraster, mask=mraster == ndv1, fill_value=ndv1)
        mmin = mraster.min()
        mraster = block_reduce(mraster, blocksize, func=how)
        araster = gdalnumeric.LoadFile(alignraster)
        araster = np.ma.masked_array(araster, mask=araster == ndv2, fill_value=ndv2)
        amin = araster.min()
        if geot1[0] <= geot2[0]:
            row3, mcol = map_pixel(geot2[0], geot2[3], geot1[1] *blocksize[0],
                                   geot1[-1]*blocksize[1], geot1[0], geot1[3])
            acol = 0
        else:
            row3, acol = map_pixel(geot1[0], geot1[3], geot2[1], geot2[-1], geot2[0], geot2[3])
            mcol = 0
        if geot1[3] <= geot2[3]:
            arow, col3 = map_pixel(geot1[0], geot1[3], geot2[1], geot2[-1], geot2[0], geot2[3])
            mrow = 0
        else:
            mrow, col3 = map_pixel(geot2[0], geot2[3], geot1[1] *blocksize[0],
                                   geot1[-1]*blocksize[1], geot1[0], geot1[3])
            arow = 0
        '''
        col3,row3 = map_pixel(geot1[0], geot1[3], geot2[1],geot2[-1], geot2[0], geot2[3])
        col3 = max(0,col3)
        row3 = max(0,row3)
        araster = araster[row3:,col3:]
        col3,row3 = map_pixel(geot2[0], geot2[3], geot1[1] *blocksize[0],
                              geot1[-1]*blocksize[1], geot1[0], geot1[3])
        col3 = max(0,abs(col3))
        row3 = max(0,np.abs(row3))
        mraster = mraster[row3:,col3:]
        '''
        mraster = mraster[mrow:, mcol:]
        araster = araster[arow:, acol:]
        if cxsize and cysize:
            araster = araster[:cysize, :cxsize]
            mraster = mraster[:cysize, :cxsize]
        else:
            rows = min(araster.shape[0], mraster.shape[0])
            cols = min(araster.shape[1], mraster.shape[1])
            araster = araster[:rows, :cols]
            mraster = mraster[:rows, :cols]
        #mraster = mraster[row3:rows+row3,col3:cols+col3]
        if masked:
            mraster = np.ma.masked_array(mraster, mask=mraster < mmin, fill_value=ndv1)
            araster = np.ma.masked_array(araster, mask=araster < amin, fill_value=ndv2)
        geot = (max(geot1[0], geot2[0]), geot1[1]*blocksize[0], geot1[2],
                min(geot1[3], geot2[3]), geot1[4], geot1[-1]*blocksize[1])
        return (mraster, araster, geot)
    else:
        print("Rasters need to be in same projection")
        return (-1, -1, -1)