code: def bilinear(arr, x, y):
      
      # for now, only 2x2 arrays
      assert arr.shape == (2, 2)
      ulv, urv, llv, lrv = arr[0:2, 0:2].flatten().tolist()
  
      # not valid if not on unit square
      assert 0.0 <= x <= 1.0
      assert 0.0 <= y <= 1.0
  
      if hasattr(arr, 'count') and arr.count() != 4:
          # a masked array with at least one nodata
          # fall back to nearest neighbor
          val = arr[int(round(1 - y)), int(round(x))]
          if val is masked:
              return None
          else:
              return asscalar(val)
  
      # bilinear interp on unit square
      return ((llv * (1 - x) * (1 - y)) +
              (lrv * x * (1 - y)) +
              (ulv * (1 - x) * y) +
              (urv * x * y))
docs: Given a 2x2 array, an x, and y, treat center points as a unit square
      return the value for the fractional row/col
      using bilinear interpolation between the cells
  
          +---+---+
          | A | B |      +----+
          +---+---+  =>  |    |
          | C | D |      +----+
          +---+---+
  
          e.g.: Center of A is at (0, 1) on unit square, D is at (1, 0), etc