code: def bilinear(px, py, band_array, gt):
      
      #import malib
      #band_array = malib.checkma(band_array)
      ndv = band_array.fill_value
      ny, nx = band_array.shape
      # Half raster cell widths
      hx = gt[1]/2.0
      hy = gt[5]/2.0
      # Calculate raster lower bound indices from point
      fx =(px -(gt[0] + hx))/gt[1]
      fy =(py -(gt[3] + hy))/gt[5]
      ix1 = int(np.floor(fx))
      iy1 = int(np.floor(fy))
      # Special case where point is on upper bounds
      if fx == float(nx - 1):
          ix1 -= 1
      if fy == float(ny - 1):
          iy1 -= 1
      # Upper bound indices on raster
      ix2 = ix1 + 1
      iy2 = iy1 + 1
      # Test array bounds to ensure point is within raster midpoints
      if(ix1 < 0) or(iy1 < 0) or(ix2 > nx - 1) or(iy2 > ny - 1):
          return ndv
      # Calculate differences from point to bounding raster midpoints
      dx1 = px -(gt[0] + ix1*gt[1] + hx)
      dy1 = py -(gt[3] + iy1*gt[5] + hy)
      dx2 =(gt[0] + ix2*gt[1] + hx) - px
      dy2 =(gt[3] + iy2*gt[5] + hy) - py
      # Use the differences to weigh the four raster values
      div = gt[1]*gt[5]
      return(band_array[iy1,ix1]*dx2*dy2/div +
              band_array[iy1,ix2]*dx1*dy2/div +
              band_array[iy2,ix1]*dx2*dy1/div +
              band_array[iy2,ix2]*dx1*dy1/div)
docs: Bilinear interpolated point at(px, py) on band_array