docs: 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)