Align two rasters so that data overlaps by geographical location
    Usage: (alignedraster_o, alignedraster_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_georasters(raster, alignraster, how=np.mean, cxsize=None, cysize=None):
    '''
    Align two rasters so that data overlaps by geographical location
    Usage: (alignedraster_o, alignedraster_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) = (raster.nodata_value, raster.shape[1],
                                                             raster.shape[0], raster.geot,
                                                             raster.projection, raster.datatype)
    (ndv2, xsize2, ysize2, geot2, projection2, datatype2) = (alignraster.nodata_value,
                                                             alignraster.shape[1],
                                                             alignraster.shape[0],
                                                             alignraster.geot,
                                                             alignraster.projection,
                                                             alignraster.datatype)
    if projection1.ExportToMICoordSys() == projection2.ExportToMICoordSys():
        blocksize = (np.round(max(geot2[1]/geot1[1], 1)).astype(np.int), np.round(max(geot2[-1]/geot1[-1], 1)).astype(np.int))
        mraster = raster.raster
        mmin = mraster.min()
        if block_reduce != (1, 1):
            mraster = block_reduce(mraster, blocksize, func=how)
        blocksize = (np.round(max(geot1[1]/geot2[1], 1)).astype(np.int), np.round(max(geot1[-1]/geot2[-1], 1)).astype(np.int))
        araster = alignraster.raster
        amin = araster.min()
        if block_reduce != (1, 1):
            araster = block_reduce(araster, blocksize, func=how)
        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
        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 = 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])
        mraster = GeoRaster(mraster, geot, projection=projection1,
                            nodata_value=ndv1, datatype=datatype1)
        araster = GeoRaster(araster, geot, projection=projection2,
                            nodata_value=ndv2, datatype=datatype2)
        return (mraster, araster)
    else:
        print("Rasters need to be in same projection")
        return (-1, -1)