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