///|
pub fn nearest_sample(
  samples : ArrayView[FieldSample],
  point : GeoPoint,
) -> FieldSample? {
  let mut result : FieldSample? = None
  let mut best = 1.0e308
  for sample in samples {
    let distance = haversine_distance_m(sample.point, point)
    if is_finite(sample.value) && distance < best {
      best = distance
      result = Some(sample)
    }
  }
  result
}

///|
pub fn k_nearest_samples(
  samples : ArrayView[FieldSample],
  point : GeoPoint,
  k : Int,
) -> Array[FieldSample] {
  let ordered = samples_sort_by_distance(samples, point)
  let output : Array[FieldSample] = []
  for i in 0.. SpatialSummary {
  let mut count = 0
  let mut total = 0.0
  let mut weighted = 0.0
  let mut weight_total = 0.0
  let mut minimum = 1.0e308
  let mut maximum = 0.0
  for sample in samples {
    let distance = haversine_distance_m(point, sample.point)
    if is_finite(sample.value) && (radius_m <= 0.0 || distance <= radius_m) {
      count += 1
      total += sample.value
      if distance < minimum {
        minimum = distance
      }
      if distance > maximum {
        maximum = distance
      }
      let weight = 1.0 / Double::max(1.0, distance)
      weighted += weight * sample.value
      weight_total += weight
    }
  }
  {
    neighbours: count,
    mean: if count == 0 {
      0.0
    } else {
      total / Double::from_int(count)
    },
    weighted_mean: if weight_total == 0.0 {
      0.0
    } else {
      weighted / weight_total
    },
    minimum_distance_m: if count == 0 {
      0.0
    } else {
      minimum
    },
    maximum_distance_m: maximum,
  }
}

///|
pub fn local_mean(
  samples : ArrayView[FieldSample],
  point : GeoPoint,
  radius_m : Double,
) -> Double {
  spatial_summary(samples, point, radius_m).mean
}

///|
pub fn local_weighted_mean(
  samples : ArrayView[FieldSample],
  point : GeoPoint,
  radius_m : Double,
) -> Double {
  spatial_summary(samples, point, radius_m).weighted_mean
}

///|
pub fn sample_density(
  samples : ArrayView[FieldSample],
  point : GeoPoint,
  radius_m : Double,
) -> Double {
  if radius_m <= 0.0 {
    0.0
  } else {
    Double::from_int(spatial_summary(samples, point, radius_m).neighbours) /
    (PI * radius_m * radius_m)
  }
}

///|
pub fn spatial_outlier_samples(
  samples : ArrayView[FieldSample],
  radius_m : Double,
  threshold : Double,
) -> Array[FieldSample] {
  let output : Array[FieldSample] = []
  for sample in samples {
    let neighbours = spatial_summary(samples, sample.point, radius_m)
    if neighbours.neighbours > 1 &&
      (sample.value - neighbours.mean).abs() > threshold {
      output.push(sample)
    }
  }
  output
}

///|
pub fn grid_sample_points(
  grid : Grid,
  points : ArrayView[GeoPoint],
) -> Array[PointEstimate] {
  let output : Array[PointEstimate] = []
  for point in points {
    match grid.sample_bilinear(point) {
      Some(value) => output.push({ point, value, neighbours: 4 })
      None => output.push({ point, value: 0.0 / 0.0, neighbours: 0 })
    }
  }
  output
}

///|
pub fn grid_sample_field(
  grid : Grid,
  samples : ArrayView[FieldSample],
) -> Array[PointEstimate] {
  let points : Array[GeoPoint] = []
  for sample in samples {
    points.push(sample.point)
  }
  grid_sample_points(grid, points)
}

///|
pub fn spatial_join_nearest(
  left : ArrayView[FieldSample],
  right : ArrayView[FieldSample],
  max_distance_m : Double,
) -> Array[(FieldSample, FieldSample)] {
  let output : Array[(FieldSample, FieldSample)] = []
  for sample in left {
    match nearest_sample(right, sample.point) {
      Some(candidate) =>
        if max_distance_m <= 0.0 ||
          haversine_distance_m(sample.point, candidate.point) <= max_distance_m {
          output.push((sample, candidate))
        }
      None => ()
    }
  }
  output
}