///|
pub fn kernel3x3(values~ : Array[Double], divisor? : Double = 1.0) -> Kernel3x3 {
  { values, divisor: if divisor == 0.0 { 1.0 } else { divisor } }
}

///|
pub fn Grid::raster_stats(self : Grid) -> RasterStats {
  let s = self.summary()
  {
    valid: s.count,
    missing: s.missing,
    minimum: s.minimum,
    maximum: s.maximum,
    mean: s.mean,
    standard_deviation: s.variance.sqrt(),
    area_m2: grid_area_m2(self),
  }
}

///|
pub fn Grid::resample(self : Grid, rows : Int, cols : Int) -> Grid {
  let safe_rows = Int::max(1, rows)
  let safe_cols = Int::max(1, cols)
  let new_lat_step = if safe_rows == 1 {
    0.0
  } else {
    self.lat_step *
    Double::from_int(self.rows - 1) /
    Double::from_int(safe_rows - 1)
  }
  let new_lon_step = if safe_cols == 1 {
    0.0
  } else {
    self.lon_step *
    Double::from_int(self.cols - 1) /
    Double::from_int(safe_cols - 1)
  }
  let values : Array[Double] = []
  for row in 0.. values.push(value)
        None => values.push(0.0 / 0.0)
      }
    }
  }
  {
    origin_lat: self.origin_lat,
    origin_lon: self.origin_lon,
    lat_step: new_lat_step,
    lon_step: new_lon_step,
    rows: safe_rows,
    cols: safe_cols,
    values,
  }
}

///|
pub fn Grid::resample_nearest(self : Grid, rows : Int, cols : Int) -> Grid {
  let safe_rows = Int::max(1, rows)
  let safe_cols = Int::max(1, cols)
  let lat_scale = if safe_rows == 1 {
    0.0
  } else {
    Double::from_int(self.rows - 1) / Double::from_int(safe_rows - 1)
  }
  let lon_scale = if safe_cols == 1 {
    0.0
  } else {
    Double::from_int(self.cols - 1) / Double::from_int(safe_cols - 1)
  }
  let values : Array[Double] = []
  for row in 0.. values.push(value)
        None => values.push(0.0 / 0.0)
      }
    }
  }
  {
    origin_lat: self.origin_lat,
    origin_lon: self.origin_lon,
    lat_step: self.lat_step * lat_scale,
    lon_step: self.lon_step * lon_scale,
    rows: safe_rows,
    cols: safe_cols,
    values,
  }
}

///|
fn kernel_value(kernel : Kernel3x3, row : Int, col : Int) -> Double {
  if row < 0 || row >= 3 || col < 0 || col >= 3 || kernel.values.length() != 9 {
    0.0
  } else {
    kernel.values[row * 3 + col] / kernel.divisor
  }
}

///|
pub fn Grid::convolve(self : Grid, kernel : Kernel3x3) -> Grid {
  let output : Array[Double] = []
  for row in 0..
              if is_finite(value) {
                total += value * kernel_value(kernel, kr, kc)
              }
            None => ()
          }
        }
      }
      output.push(total)
    }
  }
  { ..self, values: output }
}

///|
pub fn Grid::gaussian_blur(self : Grid, radius : Int) -> Grid {
  let safe = Int::max(0, Int::min(radius, 6))
  let size = safe * 2 + 1
  let weights : Array[Double] = []
  for row in 0..
              if is_finite(sample) {
                let weight = weights[kr * size + kc]
                value += sample * weight
                weight_total += weight
              }
            None => ()
          }
        }
      }
      output.push(
        if weight_total == 0.0 {
          0.0 / 0.0
        } else {
          value / weight_total
        },
      )
    }
  }
  { ..self, values: output }
}

///|
pub fn Grid::sobel_north(self : Grid) -> Grid {
  self.convolve(
    kernel3x3(values=[-1.0, -2.0, -1.0, 0.0, 0.0, 0.0, 1.0, 2.0, 1.0]),
  )
}

///|
pub fn Grid::sobel_east(self : Grid) -> Grid {
  self.convolve(
    kernel3x3(values=[-1.0, 0.0, 1.0, -2.0, 0.0, 2.0, -1.0, 0.0, 1.0]),
  )
}

///|
pub fn Grid::edge_strength(self : Grid) -> Grid {
  let north = self.sobel_north()
  let east = self.sobel_east()
  let output : Array[Double] = []
  for i in 0.. Grid {
  let output : Array[Double] = []
  for value in self.values {
    if !is_finite(value) || classes.length() == 0 {
      output.push(value)
    } else {
      let mut index = 0
      for threshold in breaks {
        if value >= threshold {
          index += 1
        }
      }
      output.push(classes[Int::min(index, classes.length() - 1)])
    }
  }
  { ..self, values: output }
}

///|
pub fn Grid::mask(
  self : Grid,
  predicate : (GridCell) -> Bool,
  replacement : Double,
) -> Grid {
  self.map_cells(fn(cell) {
    if predicate(cell) {
      cell.value
    } else {
      replacement
    }
  })
}

///|
pub fn Grid::pad(self : Grid, border : Int, value : Double) -> Grid {
  let safe = Int::max(0, border)
  let rows = self.rows + 2 * safe
  let cols = self.cols + 2 * safe
  let values : Array[Double] = []
  for row in 0.. values.push(sample)
        None => values.push(value)
      }
    }
  }
  {
    origin_lat: self.origin_lat - self.lat_step * Double::from_int(safe),
    origin_lon: self.origin_lon - self.lon_step * Double::from_int(safe),
    lat_step: self.lat_step,
    lon_step: self.lon_step,
    rows,
    cols,
    values,
  }
}

///|
pub fn Grid::crop(
  self : Grid,
  first_row : Int,
  first_col : Int,
  rows : Int,
  cols : Int,
) -> Grid {
  let safe_rows = Int::max(0, rows)
  let safe_cols = Int::max(0, cols)
  let values : Array[Double] = []
  for row in 0.. values.push(value)
        None => values.push(0.0 / 0.0)
      }
    }
  }
  {
    origin_lat: self.origin_lat + self.lat_step * Double::from_int(first_row),
    origin_lon: self.origin_lon + self.lon_step * Double::from_int(first_col),
    lat_step: self.lat_step,
    lon_step: self.lon_step,
    rows: safe_rows,
    cols: safe_cols,
    values,
  }
}

///|
pub fn Grid::valid_fraction(self : Grid) -> Double {
  if self.values.length() == 0 {
    0.0
  } else {
    Double::from_int(count_finite(self.values)) /
    Double::from_int(self.values.length())
  }
}

///|
pub fn Grid::value_range(self : Grid) -> NumericRange {
  let s = self.summary()
  { minimum: s.minimum, maximum: s.maximum }
}

///|
pub fn Grid::threshold_mask(self : Grid, range : NumericRange) -> Grid {
  self.map(fn(value) { if range_contains(range, value) { 1.0 } else { 0.0 } })
}

///|
pub fn Grid::quantize(self : Grid, step : Double) -> Grid {
  if step <= 0.0 {
    self
  } else {
    self.map(fn(value) {
      if is_finite(value) {
        (value / step).round() * step
      } else {
        value
      }
    })
  }
}

///|
pub fn Grid::connected_components(
  self : Grid,
  predicate : (Double) -> Bool,
) -> Array[RegionSummary] {
  let labels : Array[Int] = []
  for _ in self.values {
    labels.push(0)
  }
  let regions : Array[RegionSummary] = []
  let mut label = 0
  for row in 0.. maximum {
              maximum = value
            }
            total += value
            cells += 1
          }
          let neighbours = [
            current - self.cols,
            current + self.cols,
            current - 1,
            current + 1,
          ]
          for next in neighbours {
            if next >= 0 && next < self.values.length() {
              let nr = next / self.cols
              let nc = next % self.cols
              if (nr == cr || nc == cc) &&
                labels[next] == 0 &&
                predicate(self.values[next]) {
                labels[next] = label
                queue.push(next)
              }
            }
          }
        }
        regions.push({
          label,
          cells,
          mean: if cells == 0 {
            0.0
          } else {
            total / Double::from_int(cells)
          },
          minimum,
          maximum,
          area_m2: Double::from_int(cells) * grid_cell_area_m2(self, row),
        })
      }
    }
  }
  regions
}

///|
pub fn idw_estimate(
  samples : ArrayView[FieldSample],
  point : GeoPoint,
  power? : Double = 2.0,
  radius_m? : Double = 0.0,
) -> PointEstimate {
  let mut weighted = 0.0
  let mut total_weight = 0.0
  let mut neighbours = 0
  for sample in samples {
    if is_finite(sample.value) {
      let distance = haversine_distance_m(sample.point, point)
      if distance == 0.0 {
        return { point, value: sample.value, neighbours: 1 }
      }
      if radius_m == 0.0 || distance <= radius_m {
        let weight = 1.0 / @math.pow(distance, power)
        weighted += weight * sample.value
        total_weight += weight
        neighbours += 1
      }
    }
  }
  {
    point,
    value: if total_weight == 0.0 {
      0.0 / 0.0
    } else {
      weighted / total_weight
    },
    neighbours,
  }
}

///|
pub fn grid_from_samples_idw(
  grid : Grid,
  samples : ArrayView[FieldSample],
  power? : Double = 2.0,
  radius_m? : Double = 0.0,
) -> Grid {
  let output : Array[Double] = []
  for row in 0..
          output.push(idw_estimate(samples, point, power~, radius_m~).value)
        None => output.push(0.0 / 0.0)
      }
    }
  }
  { ..grid, values: output }
}

///|
pub fn grid_to_samples(
  grid : Grid,
  name? : String = "grid",
) -> Array[FieldSample] {
  let output : Array[FieldSample] = []
  for cell in grid.cells() {
    output.push({ point: cell.point, value: cell.value, name })
  }
  output
}

///|
pub fn Grid::grid_region_summaries(
  self : Grid,
  labels : Grid,
) -> Array[RegionSummary] {
  let groups : Map[Int, Array[Double]] = Map([])
  for i in 0.. existing
        None => []
      }
      values.push(self.values[i])
      groups[label] = values
    }
  }
  let output : Array[RegionSummary] = []
  for label, values in groups {
    let s = summary(values)
    output.push({
      label,
      cells: s.count,
      mean: s.mean,
      minimum: s.minimum,
      maximum: s.maximum,
      area_m2: Double::from_int(s.count) * grid_cell_area_m2(self, 0),
    })
  }
  output
}