diff --git a/clayground.py b/clayground.py
deleted file mode 100644
index a73fec94..00000000
--- a/clayground.py
+++ /dev/null
@@ -1,132 +0,0 @@
-import lancedb
-import matplotlib.pyplot as plt
-import rasterio as rio
-import streamlit as st
-from rasterio.plot import show
-
-st.set_page_config(layout="wide")
-
-
-# Get preferrred chips
-def get_unique_chips(tbl):
-    chips = [
-        {"tile": "17MNP", "idx": "0271", "year": 2023},
-        {"tile": "19HGU", "idx": "0033", "year": 2018},
-        {"tile": "33NVB", "idx": "0393", "year": 2020},
-        {"tile": "21JVJ", "idx": "0100", "year": 2020},
-        {"tile": "34KHD", "idx": "0080", "year": 2018},
-        {"tile": "19JCF", "idx": "0215", "year": 2023},
-        {"tile": "20HMK", "idx": "0100", "year": 2020},
-        {"tile": "37MFT", "idx": "0313", "year": 2023},
-        {"tile": "49KHR", "idx": "0020", "year": 2017},
-        {"tile": "55LBC", "idx": "0075", "year": 2022},
-    ]
-
-    tile_filter = " OR ".join(
-        [
-            f"(tile == '{chip['tile']}' "
-            f"AND idx == '{chip['idx']}') "
-            f"AND year == {chip['year']}"
-            for chip in chips
-        ]
-    )
-    result = tbl.search().where(tile_filter, prefilter=True).to_pandas()
-    return result
-
-
-# Load embeddings
-@st.cache_resource()
-def connect_to_database():
-    db = lancedb.connect("nbs/embeddings")
-    tbl = db.open_table("clay-v001")
-    return tbl
-
-
-@st.cache_resource()
-def show_samples(_tbl):
-    df = get_unique_chips(_tbl)
-    # df = _tbl.head(10).to_pandas()
-    # sample 100 random rows
-    # samples = df.sample(100).to_dict("records")
-    samples = df.to_dict("records")
-
-    cols = st.columns(10)
-    options = {}
-    for idx, sample in enumerate(samples):
-        path = sample["path"]
-        rgb_chip = rio.open(path).read(indexes=[3, 2, 1])
-        rgb_chip = (rgb_chip - rgb_chip.min()) / (rgb_chip.max() - rgb_chip.min())
-        with cols[idx % 10]:
-            st.caption(f"{sample['tile']}-{sample['date']}-{sample['idx']}")
-            show(rgb_chip)
-            plt.axis("off")
-            st.pyplot(plt)
-
-        options[f"{sample['tile']}-{sample['idx']}"] = {
-            "vector": sample["vector"],
-            "tile": sample["tile"],
-            "year": sample["year"],
-        }
-
-    return options
-
-
-# Function to find similar vectors
-@st.cache_data()
-def find_similar_vectors(_tbl, query):
-    # tile, year = query["tile"], query["year"]
-    # filter = f"tile != '{tile}'"
-    result = (
-        _tbl.search(query=query["vector"], vector_column_name="vector")
-        .metric("cosine")
-        # .where(filter, prefilter=True)
-        .limit(10)
-        .to_pandas()
-    )
-    # st.dataframe(result)
-    cols = st.columns(10)
-    for idx, row in result.iterrows():
-        path = row["path"]
-        rgb_chip = rio.open(path).read(indexes=[3, 2, 1])
-        rgb_chip = (rgb_chip - rgb_chip.min()) / (rgb_chip.max() - rgb_chip.min())
-        with cols[idx % 10]:
-            st.caption(f"{row['tile']}-{row['date']}-{row['idx']}")
-            show(rgb_chip)
-            plt.axis("off")
-            st.pyplot(plt)
-
-
-# Main app
-def main():
-    st.title("Clayground")
-
-    tbl = connect_to_database()
-    options = show_samples(tbl)
-
-    # UI to select an embedding
-    with st.sidebar:
-        selection = st.selectbox("Select a chip", options=options.keys())
-
-        arithmetic = st.toggle("Arithmetic", False)
-        if arithmetic:
-            multiselect = st.multiselect(
-                "Select multiple chips", options=options.keys(), default=[]
-            )
-
-        submit = st.button("Submit")
-
-    if submit and not arithmetic:
-        query = options[selection]
-        find_similar_vectors(tbl, query)
-
-    if submit and arithmetic and len(multiselect) > 1:
-        st.write("Selected:", multiselect)
-        v1 = options[multiselect[0]]
-        v2 = options[multiselect[1]]
-        v3 = (v1["vector"] + v2["vector"]) / 2
-
-        find_similar_vectors(tbl, {"vector": v3})
-
-
-if __name__ == "__main__":
-    main()
diff --git a/scripts/pipeline/batch/Dockerfile b/scripts/pipeline/batch/Dockerfile
deleted file mode 100644
index 24e4a755..00000000
--- a/scripts/pipeline/batch/Dockerfile
+++ /dev/null
@@ -1,8 +0,0 @@
-# List of aws images https://github.com/aws/deep-learning-containers/blob/master/available_images.md
-FROM mcr.microsoft.com/planetary-computer/python
-
-# Fetch-and-run setup.
-# https://aws.amazon.com/blogs/compute/creating-a-simple-fetch-and-run-aws-batch-job/
-ADD fetch_and_run.sh /usr/local/bin/fetch_and_run.sh
-WORKDIR /tmp
-ENTRYPOINT ["/usr/local/bin/fetch_and_run.sh"]
diff --git a/scripts/pipeline/batch/compute-environment.json b/scripts/pipeline/batch/compute-environment.json
deleted file mode 100644
index 4e40e4c0..00000000
--- a/scripts/pipeline/batch/compute-environment.json
+++ /dev/null
@@ -1,42 +0,0 @@
-{
-  "computeEnvironmentName": "fetch-and-run",
-  "computeEnvironmentArn": "arn:aws:batch:us-east-1:939659358055:compute-environment/fetch-and-run",
-  "ecsClusterArn": "arn:aws:ecs:us-east-1:939659358055:cluster/AWSBatch-fetch-and-run-6c59b338-7c7a-3d17-9cae-48fa8cb01fd1",
-  "tags": {},
-  "type": "MANAGED",
-  "state": "ENABLED",
-  "status": "VALID",
-  "statusReason": "ComputeEnvironment Healthy",
-  "computeResources": {
-    "type": "SPOT",
-    "allocationStrategy": "SPOT_PRICE_CAPACITY_OPTIMIZED",
-    "minvCpus": 0,
-    "maxvCpus": 256,
-    "desiredvCpus": 0,
-    "instanceTypes": [
-      "optimal"
-    ],
-    "subnets": [
-      "subnet-065db57a9bbe9db39",
-      "subnet-071e09a3c54797241",
-      "subnet-02b3216f099a85184",
-      "subnet-082a6aa9a80b134ce",
-      "subnet-021cbf7ad41e521ff",
-      "subnet-0ea444ed7a76a64ee"
-    ],
-    "securityGroupIds": [
-      "sg-076376c48cdb15a88"
-    ],
-    "instanceRole": "arn:aws:iam::939659358055:instance-profile/ecsInstanceRole",
-    "tags": {},
-    "bidPercentage": 100,
-    "ec2Configuration": [
-      {
-        "imageType": "ECS_AL2"
-      }
-    ]
-  },
-  "serviceRole": "arn:aws:iam::939659358055:role/aws-service-role/batch.amazonaws.com/AWSServiceRoleForBatch",
-  "containerOrchestrationType": "ECS",
-  "uuid": "b105b281-f736-36de-af70-450cafda628a"
-}
diff --git a/scripts/pipeline/batch/fetch_and_run.sh b/scripts/pipeline/batch/fetch_and_run.sh
deleted file mode 100755
index c34e5bde..00000000
--- a/scripts/pipeline/batch/fetch_and_run.sh
+++ /dev/null
@@ -1,128 +0,0 @@
-#!/bin/bash
-
-# Copyright 2013-2017 Amazon.com, Inc. or its affiliates. All Rights Reserved.
-#
-# Licensed under the Apache License, Version 2.0 (the "License"). You may not use this file except in compliance with the
-# License. A copy of the License is located at
-#
-# http://aws.amazon.com/apache2.0/
-#
-# or in the "LICENSE.txt" file accompanying this file. This file is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES
-# OR CONDITIONS OF ANY KIND, express or implied. See the License for the specific language governing permissions and
-# limitations under the License.
-
-# This script can help you download and run a script from S3 using aws-cli.
-# It can also download a zip file from S3 and run a script from inside.
-# See below for usage instructions.
-
-# https://raw.githubusercontent.com/awslabs/aws-batch-helpers/master/fetch-and-run/fetch_and_run.sh
-
-PATH="/srv/conda/envs/notebook/bin:/opt/conda/bin:/bin:/usr/bin:/sbin:/usr/sbin:/usr/local/bin:/usr/local/sbin"
-BASENAME="${0##*/}"
-
-usage () {
-  if [ "${#@}" -ne 0 ]; then
-    echo "* ${*}"
-    echo
-  fi
-  cat <<ENDUSAGE
-Usage:
-
-export BATCH_FILE_TYPE="script"
-export BATCH_FILE_S3_URL="s3://my-bucket/my-script"
-${BASENAME} script-from-s3 [ <script arguments> ]
-
-  - or -
-
-export BATCH_FILE_TYPE="zip"
-export BATCH_FILE_S3_URL="s3://my-bucket/my-zip"
-${BASENAME} script-from-zip [ <script arguments> ]
-ENDUSAGE
-
-  exit 2
-}
-
-# Standard function to print an error and exit with a failing return code
-error_exit () {
-  echo "${BASENAME} - ${1}" >&2
-  exit 1
-}
-
-# Check what environment variables are set
-if [ -z "${BATCH_FILE_TYPE}" ]; then
-  usage "BATCH_FILE_TYPE not set, unable to determine type (zip/script) of URL ${BATCH_FILE_S3_URL}"
-fi
-
-if [ -z "${BATCH_FILE_S3_URL}" ]; then
-  usage "BATCH_FILE_S3_URL not set. No object to download."
-fi
-
-scheme="$(echo "${BATCH_FILE_S3_URL}" | cut -d: -f1)"
-if [ "${scheme}" != "s3" ]; then
-  usage "BATCH_FILE_S3_URL must be for an S3 object; expecting URL starting with s3://"
-fi
-
-# Check that necessary programs are available
-which aws >/dev/null 2>&1 || error_exit "Unable to find AWS CLI executable."
-which unzip >/dev/null 2>&1 || error_exit "Unable to find unzip executable."
-
-# Create a temporary directory to hold the downloaded contents, and make sure
-# it's removed later, unless the user set KEEP_BATCH_FILE_CONTENTS.
-cleanup () {
-   if [ -z "${KEEP_BATCH_FILE_CONTENTS}" ] \
-     && [ -n "${TMPDIR}" ] \
-     && [ "${TMPDIR}" != "/" ]; then
-      rm -r "${TMPDIR}"
-   fi
-}
-trap 'cleanup' EXIT HUP INT QUIT TERM
-# mktemp arguments are not very portable.  We make a temporary directory with
-# portable arguments, then use a consistent filename within.
-TMPDIR="$(mktemp -d -t tmp.XXXXXXXXX)" || error_exit "Failed to create temp directory."
-TMPFILE="${TMPDIR}/batch-file-temp"
-install -m 0600 /dev/null "${TMPFILE}" || error_exit "Failed to create temp file."
-
-# Fetch and run a script
-fetch_and_run_script () {
-  # Create a temporary file and download the script
-  aws s3 cp "${BATCH_FILE_S3_URL}" - > "${TMPFILE}" || error_exit "Failed to download S3 script."
-
-  # Make the temporary file executable and run it with any given arguments
-  local script="./${1}"; shift
-  chmod u+x "${TMPFILE}" || error_exit "Failed to chmod script."
-  exec ${TMPFILE} "${@}" || error_exit "Failed to execute script."
-}
-
-# Download a zip and run a specified script from inside
-fetch_and_run_zip () {
-  # Create a temporary file and download the zip file
-  aws s3 cp "${BATCH_FILE_S3_URL}" - > "${TMPFILE}" || error_exit "Failed to download S3 zip file from ${BATCH_FILE_S3_URL}"
-
-  # Create a temporary directory and unpack the zip file
-  cd "${TMPDIR}" || error_exit "Unable to cd to temporary directory."
-  unzip -q "${TMPFILE}" || error_exit "Failed to unpack zip file."
-
-  # Use first argument as script name and pass the rest to the script
-  local script="./${1}"; shift
-  [ -r "${script}" ] || error_exit "Did not find specified script '${script}' in zip from ${BATCH_FILE_S3_URL}"
-  chmod u+x "${script}" || error_exit "Failed to chmod script."
-  exec "${script}" "${@}" || error_exit " Failed to execute script."
-}
-
-# Main - dispatch user request to appropriate function
-case ${BATCH_FILE_TYPE} in
-  zip)
-    if [ ${#@} -eq 0 ]; then
-      usage "zip format requires at least one argument - the script to run from inside"
-    fi
-    fetch_and_run_zip "${@}"
-    ;;
-
-  script)
-    fetch_and_run_script "${@}"
-    ;;
-
-  *)
-    usage "Unsupported value for BATCH_FILE_TYPE. Expected (zip/script)."
-    ;;
-esac
diff --git a/scripts/pipeline/batch/job-definition.json b/scripts/pipeline/batch/job-definition.json
deleted file mode 100644
index f1cc0504..00000000
--- a/scripts/pipeline/batch/job-definition.json
+++ /dev/null
@@ -1,59 +0,0 @@
-{
-  "jobDefinitionName": "fetch-and-run",
-  "jobDefinitionArn": "arn:aws:batch:us-east-1:939659358055:job-definition/fetch-and-run:1",
-  "revision": 1,
-  "status": "ACTIVE",
-  "type": "container",
-  "parameters": {},
-  "retryStrategy": {
-    "attempts": 5,
-    "evaluateOnExit": [
-      {
-        "onStatusReason": "Host EC2*",
-        "action": "retry"
-      },
-      {
-        "onReason": "*",
-        "action": "exit"
-      }
-    ]
-  },
-  "containerProperties": {
-    "image": "939659358055.dkr.ecr.us-east-1.amazonaws.com/fetch-and-run:latest",
-    "command": [
-      "datacube.py",
-      "--help"
-    ],
-    "jobRoleArn": "arn:aws:iam::939659358055:role/ecsS3FullAccess",
-    "executionRoleArn": "arn:aws:iam::939659358055:role/ecsTaskExecutionRole",
-    "volumes": [],
-    "environment": [],
-    "mountPoints": [],
-    "ulimits": [],
-    "resourceRequirements": [
-      {
-        "value": "6",
-        "type": "VCPU"
-      },
-      {
-        "value": "12000",
-        "type": "MEMORY"
-      }
-    ],
-    "linuxParameters": {
-      "devices": [],
-      "tmpfs": []
-    },
-    "logConfiguration": {
-      "logDriver": "awslogs",
-      "options": {},
-      "secretOptions": []
-    },
-    "secrets": []
-  },
-  "tags": {},
-  "platformCapabilities": [
-    "EC2"
-  ],
-  "containerOrchestrationType": "ECS"
-}
diff --git a/scripts/pipeline/batch/submit.py b/scripts/pipeline/batch/submit.py
deleted file mode 100644
index cbc568ab..00000000
--- a/scripts/pipeline/batch/submit.py
+++ /dev/null
@@ -1,40 +0,0 @@
-import os
-
-import boto3
-
-batch = boto3.client("batch", region_name="us-east-1")
-
-NR_OF_TILES_IN_SAMPLE_FILE = 2728
-
-PC_KEY = os.environ["PC_SDK_SUBSCRIPTION_KEY"]
-
-
-job = {
-    "jobName": "fetch-and-run-datacube",
-    "jobQueue": "fetch-and-run",
-    "jobDefinition": "fetch-and-run",
-    "containerOverrides": {
-        "command": [
-            "datacube.py",
-            "--bucket",
-            "clay-tiles-04-sample-v02",
-            "--sample",
-            "https://clay-mgrs-samples.s3.amazonaws.com/mgrs_sample_v02.fgb",
-        ],
-        "environment": [
-            {"name": "BATCH_FILE_TYPE", "value": "zip"},
-            {
-                "name": "BATCH_FILE_S3_URL",
-                "value": "s3://clay-fetch-and-run-packages/batch-fetch-and-run.zip",
-            },
-            {"name": "PC_SDK_SUBSCRIPTION_KEY", "value": f"{PC_KEY}"},
-        ],
-        "resourceRequirements": [
-            {"type": "MEMORY", "value": "15500"},
-            {"type": "VCPU", "value": "4"},
-        ],
-    },
-    "arrayProperties": {"size": NR_OF_TILES_IN_SAMPLE_FILE},
-}
-
-print(batch.submit_job(**job))
diff --git a/scripts/pipeline/datacube.py b/scripts/pipeline/datacube.py
deleted file mode 100755
index 43913eba..00000000
--- a/scripts/pipeline/datacube.py
+++ /dev/null
@@ -1,521 +0,0 @@
-#!/usr/bin/env python3
-"""
-STAC Data Processing Script
-
-This Python script processes Sentinel-2, Sentinel-1, and Copernicus DEM
-(Digital Elevation Model) data. It utilizes Microsoft's Planetary Computer API
-for data retrieval and manipulation.
-
-Constants:
-- STAC_API: Planetary Computer API endpoint
-- S2_BANDS: Bands used in Sentinel-2 data processing
-
-Functions:
-- get_surrounding_days(reference, interval_days):
-      Get the week range for a given date.
-- search_sentinel2(date_range, aoi, cloud_cover_percentage, nodata_pixel_percentage):
-      Search for Sentinel-2 items within a given date range and area of interest.
-- search_sentinel1(bbox, catalog, date_range):
-      Search for Sentinel-1 items within a given bounding box, STAC catalog,
-      and date range.
-- search_dem(bbox, catalog):
-      Search for DEM items within a given bounding box.
-- make_datasets(s2_item, s1_items, dem_items, resolution):
-      Create xarray Datasets for Sentinel-2, Sentinel-1, and DEM data.
-- process(aoi, year, resolution, cloud_cover_percentage, nodata_pixel_percentage):
-      Process Sentinel-2, Sentinel-1, and DEM data for a specified time range,
-      area of interest, and resolution.
-"""
-
-import logging
-import os
-import random
-from datetime import timedelta
-
-import click
-import geopandas as gpd
-import numpy as np
-import planetary_computer as pc
-import pystac_client
-import stackstac
-import xarray as xr
-from pystac import ItemCollection
-from shapely.geometry import box
-from tile import tiler
-
-STAC_API = "https://planetarycomputer.microsoft.com/api/stac/v1"
-S2_BANDS = ["B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12", "SCL"]
-SPATIAL_RESOLUTION = 10
-CLOUD_COVER_PERCENTAGE = 50
-NODATA_PIXEL_PERCENTAGE = 20
-NODATA = 0
-S1_MATCH_ATTEMPTS = 40
-DATES_PER_LOCATION = 4
-
-logger = logging.getLogger("datacube")
-hdr = logging.StreamHandler()
-formatter = logging.Formatter("[%(asctime)s] %(levelname)s: %(message)s")
-hdr.setFormatter(formatter)
-logger.addHandler(hdr)
-logger.setLevel(logging.INFO)
-
-
-def get_surrounding_days(reference, interval_days):
-    """
-    Get the days surrounding the input date.
-
-    Parameters:
-    - reference (datetime): The reference datetime.
-    - interval_days (int): The number of days to search ahead and back
-
-    Returns:
-    - str: A string representing the start and end dates of the date interval in the
-        format 'start_date/end_date'.
-    """
-    start = reference - timedelta(days=interval_days)
-    end = reference + timedelta(days=interval_days)
-    return f"{start.date()}/{end.date()}"
-
-
-def search_sentinel2(
-    catalog, date_range, aoi, cloud_cover_percentage, nodata_pixel_percentage, index=0
-):
-    """
-    Search for Sentinel-2 items within a given date range and area of interest (AOI)
-    with specified conditions.
-
-    Parameters:
-    - catalog (pystac.Catalog): STAC catalog containing Sentinel-2 items.
-    - date_range (str): The date range in the format 'start_date/end_date'.
-    - aoi (shapely.geometry.base.BaseGeometry): Geometry object for an Area of
-        Interest (AOI).
-    - cloud_cover_percentage (int): Maximum acceptable cloud cover percentage
-        for Sentinel-2 images.
-    - nodata_pixel_percentage (int): Maximum acceptable percentage of nodata
-        pixels in Sentinel-2 images.
-    - index: Which of the found scenes to select
-
-    Returns:
-    - tuple: A tuple containing the STAC catalog, Sentinel-2 items, and the
-        bounding box (bbox)
-
-    Note:
-    The function filters Sentinel-2 items based on the specified conditions
-    such as geometry, date, cloud cover, and nodata pixel percentage. Only one
-    result with the least cloud cover will be returned. The result is returned
-    as a tuple containing the STAC catalog, Sentinel-2 items, the bounding box
-    of the first item, and an EPSG code for the coordinate reference system.
-    """
-    search: pystac_client.item_search.ItemSearch = catalog.search(
-        filter_lang="cql2-json",
-        filter={
-            "op": "and",
-            "args": [
-                {
-                    "op": "s_intersects",
-                    "args": [{"property": "geometry"}, aoi.centroid.__geo_interface__],
-                },
-                {"op": "anyinteracts", "args": [{"property": "datetime"}, date_range]},
-                {"op": "=", "args": [{"property": "collection"}, "sentinel-2-l2a"]},
-                {
-                    "op": "<=",
-                    "args": [{"property": "eo:cloud_cover"}, cloud_cover_percentage],
-                },
-                {
-                    "op": "<=",
-                    "args": [
-                        {"property": "s2:nodata_pixel_percentage"},
-                        nodata_pixel_percentage,
-                    ],
-                },
-            ],
-        },
-    )
-
-    s2_items = search.item_collection()
-    logger.info(f"Found {len(s2_items)} Sentinel-2 items")
-    if not len(s2_items):
-        return None, None
-
-    s2_items_gdf = gpd.GeoDataFrame.from_features(s2_items.to_dict())
-
-    least_clouds = s2_items_gdf.sort_values(by=["eo:cloud_cover"], ascending=True).iloc[
-        index
-    ]
-
-    # Get the datetime for the filtered Sentinel 2 dataframe
-    # containing the least nodata and least cloudy scene
-    for item in s2_items:
-        if item.properties["datetime"] == least_clouds["datetime"]:
-            s2_item = item
-            break
-
-    bbox = least_clouds.geometry.bounds
-
-    epsg = s2_item.properties["proj:epsg"]
-    logger.info(f"EPSG code based on Sentinel-2 item: {epsg}")
-
-    return s2_item, bbox
-
-
-def search_sentinel1(bbox, catalog, date_range):
-    """
-    Search for Sentinel-1 items within a given bounding box (bbox), STAC
-    catalog, and date range.
-
-    Parameters:
-    - bbox (tuple): Bounding box coordinates in the format
-        (minx, miny, maxx, maxy).
-    - catalog (pystac.Catalog): STAC catalog containing Sentinel-1 items.
-    - date_range (str): The date range in the format 'start_date/end_date'.
-
-    Returns:
-    - pystac.Collection: A collection of Sentinel-1 items filtered by specified
-        conditions.
-
-    Note:
-    This function retrieves Sentinel-1 items from the catalog that intersect
-    with the given bounding box and fall within the provided time window. The
-    function filters items based on orbit state and returns the collection of
-    Sentinel-1 items that meet the defined criteria.
-    """
-    # Create poly geom object from the bbox
-    geom_bbox = box(*bbox)
-
-    search: pystac_client.item_search.ItemSearch = catalog.search(
-        filter_lang="cql2-json",
-        filter={
-            "op": "and",
-            "args": [
-                {
-                    "op": "s_intersects",
-                    "args": [{"property": "geometry"}, geom_bbox.__geo_interface__],
-                },
-                {"op": "anyinteracts", "args": [{"property": "datetime"}, date_range]},
-                {"op": "=", "args": [{"property": "collection"}, "sentinel-1-rtc"]},
-            ],
-        },
-    )
-    s1_items = search.item_collection()
-    logger.info(f"Found {len(s1_items)} Sentinel-1 items")
-
-    if not len(s1_items):
-        return
-    else:
-        # Add id as property to persist in gdf
-        for item in s1_items:
-            item.properties["id"] = item.id
-
-        # Try to find enough scenes with same orbit that fully overlap
-        # the S2 bbox.
-        s1_gdf = gpd.GeoDataFrame.from_features(s1_items)
-        s1_gdf["overlap"] = s1_gdf.intersection(box(*bbox)).area
-        s1_gdf = s1_gdf.sort_values(by="overlap", ascending=False)
-
-        most_overlap_orbit = s1_gdf.iloc[0]["sat:orbit_state"]
-        logger.info(f"Most overlapped orbit: {most_overlap_orbit}")
-        selected_item_ids = []
-        intersection = None
-        orbit = None
-        for index, row in s1_gdf.iterrows():
-            orbit = row["sat:orbit_state"]
-            if intersection is None and orbit == most_overlap_orbit:
-                intersection = row.geometry
-                selected_item_ids.append(row.id)
-                intersection = intersection.intersection(row.geometry)
-            elif orbit == most_overlap_orbit and not intersection.covers(geom_bbox):
-                intersection = row.geometry
-                selected_item_ids.append(row.id)
-                intersection = intersection.intersection(row.geometry)
-            elif orbit == most_overlap_orbit and intersection.covers(geom_bbox):
-                # Stop adding scenes when the bbox is fully covered.
-                break
-            else:
-                pass
-
-        s1_items = ItemCollection(
-            [item for item in s1_items if item.id in selected_item_ids]
-        )
-
-        return s1_items
-
-
-def search_dem(bbox, catalog):
-    """
-    Search for Copernicus Digital Elevation Model (DEM) items within a given
-    bounding box (bbox), STAC catalog, and Sentinel-2 items.
-
-    Parameters:
-    - bbox (tuple): Bounding box coordinates in the format
-        (minx, miny, maxx, maxy).
-    - catalog (pystac.Catalog): STAC catalog containing DEM items.
-
-    Returns:
-    - pystac.Collection: A collection of Digital Elevation Model (DEM) items
-        filtered by specified conditions.
-    """
-    search = catalog.search(collections=["cop-dem-glo-30"], bbox=bbox)
-    dem_items = search.item_collection()
-    logger.info(f"Found {len(dem_items)} DEM items")
-
-    return dem_items
-
-
-def make_datasets(s2_items, s1_items, dem_items, resolution):
-    """
-    Create xarray Datasets for Sentinel-2, Sentinel-1, and Copernicus DEM
-    data.
-
-    Parameters:
-    - s2_items (list): List of Sentinel-2 items.
-    - s1_items (list): List of Sentinel-1 items.
-    - dem_items (list): List of DEM items.
-    - resolution (int): Spatial resolution.
-
-    Returns:
-    - tuple: A tuple containing xarray Datasets for Sentinel-2, Sentinel-1,
-        and Copernicus DEM.
-    """
-    da_sen2: xr.DataArray = stackstac.stack(
-        items=s2_items,
-        assets=S2_BANDS,
-        resolution=resolution,
-        dtype=np.uint16,
-        fill_value=0,
-    )
-
-    da_sen1: xr.DataArray = stackstac.stack(
-        items=s1_items,
-        assets=["vh", "vv"],
-        epsg=int(da_sen2.epsg),
-        bounds=da_sen2.spec.bounds,
-        resolution=resolution,
-        dtype=np.float32,
-        fill_value=np.nan,
-    )
-
-    da_dem: xr.DataArray = stackstac.stack(
-        items=dem_items,
-        epsg=int(da_sen2.epsg),
-        bounds=da_sen2.spec.bounds,
-        resolution=resolution,
-        dtype=np.float32,
-        fill_value=np.nan,
-    )
-
-    da_sen1: xr.DataArray = stackstac.mosaic(da_sen1, dim="time")
-
-    da_sen1 = da_sen1.drop_vars(
-        [var for var in da_sen1.coords if var not in da_sen1.dims]
-    )
-
-    da_sen2 = da_sen2.drop_vars(
-        [var for var in da_sen2.coords if var not in da_sen2.dims]
-    ).squeeze()
-
-    del da_sen2.coords["time"]
-
-    da_dem: xr.DataArray = stackstac.mosaic(da_dem, dim="time").assign_coords(
-        {"band": ["dem"]}
-    )
-
-    da_dem = da_dem.drop_vars([var for var in da_dem.coords if var not in da_dem.dims])
-
-    return [da_sen2, da_sen1, da_dem]
-
-
-def process(
-    aoi,
-    date_range,
-    resolution,
-    cloud_cover_percentage,
-    nodata_pixel_percentage,
-):
-    """
-    Process Sentinel-2, Sentinel-1, and Copernicus DEM data for a specified
-    date_range, area of interest (AOI), resolution, EPSG code, cloud cover
-    percentage, and nodata pixel percentage.
-
-    Parameters:
-    - aoi (shapely.geometry.base.BaseGeometry): Geometry object for an Area of
-        Interest (AOI).
-    - date_range (str): Date range string to pass to the catalog search.
-    - resolution (int): Spatial resolution.
-    - cloud_cover_percentage (int): Maximum acceptable cloud cover percentage
-        for Sentinel-2 images.
-    - nodata_pixel_percentage (int): Maximum acceptable percentage of nodata
-        pixels in Sentinel-2 images.
-
-    Returns:
-    - xr.Dataset: Merged xarray Dataset containing processed data.
-    """
-    catalog = pystac_client.Client.open(STAC_API, modifier=pc.sign_inplace)
-
-    for i in range(S1_MATCH_ATTEMPTS):
-        s2_item, bbox = search_sentinel2(
-            catalog,
-            date_range,
-            aoi,
-            cloud_cover_percentage,
-            nodata_pixel_percentage,
-            index=i,
-        )
-        if not s2_item:
-            continue
-
-        surrounding_days = get_surrounding_days(s2_item.datetime, interval_days=3)
-        logger.info(f"Searching S1 in date range {surrounding_days}")
-
-        s1_items = search_sentinel1(bbox, catalog, surrounding_days)
-
-        if s1_items:
-            break
-
-    if i == S1_MATCH_ATTEMPTS - 1:
-        logger.info(
-            "No match for S1 scenes found for date range "
-            f"{date_range} after {S1_MATCH_ATTEMPTS} attempts."
-        )
-        return None, None
-
-    dem_items = search_dem(bbox, catalog)
-
-    date = s2_item.properties["datetime"][:10]
-
-    result = make_datasets(
-        s2_item,
-        s1_items,
-        dem_items,
-        resolution,
-    )
-
-    if 0 in (dat.shape[0] for dat in result):
-        logger.info("S2/S1 pixel coverages do not overlap although bounds do")
-        return None, None
-
-    return date, result
-
-
-def convert_attrs_and_coords_objects_to_str(data):
-    """
-    Convert attributes and coordinates that are objects to
-    strings.
-
-    This is required for storing the xarray in netcdf.
-    """
-    for key, coord in data.coords.items():
-        if coord.dtype == "object":
-            data.coords[key] = str(coord.values)
-
-    for key, attr in data.attrs.items():
-        data.attrs[key] = str(attr)
-
-    for key, var in data.variables.items():
-        var.attrs = {}
-
-
-@click.command()
-@click.option(
-    "--sample",
-    required=False,
-    default="https://clay-mgrs-samples.s3.amazonaws.com/mgrs_sample_v02.fgb",
-    help="Location of MGRS tile sample",
-)
-@click.option(
-    "--index",
-    required=False,
-    default=None,
-    help="Index of MGRS tile from sample file that should be processed",
-)
-@click.option(
-    "--bucket",
-    required=False,
-    default="",
-    help="Specify the bucket for where to write the data.",
-)
-@click.option(
-    "--subset",
-    required=False,
-    default="",
-    help="For debugging, subset the MGRS tile data to this pixel window."
-    "Expects a comma separated string of 4 integers.",
-    type=str,
-)
-@click.option(
-    "--localpath",
-    required=False,
-    default=None,
-    help="If specified, this path will be used to write the tiles locally"
-    "Otherwise a temp dir will be used.",
-)
-@click.option(
-    "--dateranges",
-    required=False,
-    default="",
-    type=str,
-    help="Comma separated list of date ranges, each provided as YYYY-MM-DD/YYYY-MM-DD.",
-)
-@click.option("-v", "--verbose", is_flag=True)
-def main(sample, index, subset, bucket, localpath, dateranges, verbose):  # noqa: PLR0913
-    if verbose:
-        logger.setLevel(logging.DEBUG)
-
-    if index is None:
-        index = int(os.environ.get("AWS_BATCH_JOB_ARRAY_INDEX", 0))
-    else:
-        index = int(index)
-    tiles = gpd.read_file(sample)
-    tile = tiles.iloc[index]
-    mgrs = tile["name"]
-
-    logger.info(f"Starting algorithm for MGRS tile {tile['name']} with index {index}")
-
-    if subset:
-        subset = [int(dat) for dat in subset.split(",")]
-
-    if dateranges:
-        date_ranges = dateranges.split(",")
-    else:
-        # Shuffle years, use index as seed for reproducibility but no
-        # to have the same shuffle every time.
-        date_ranges = [
-            f"{year}-01-01/{year}-12-31"
-            for year in (2017, 2018, 2019, 2020, 2021, 2022, 2023)
-        ]
-        random.seed(index)
-        random.shuffle(date_ranges)
-
-    match_count = 0
-    for date_range in date_ranges:
-        logger.info(f"Processing data for date range {date_range}")
-        date, pixels = process(
-            tile.geometry,
-            date_range,
-            SPATIAL_RESOLUTION,
-            CLOUD_COVER_PERCENTAGE,
-            NODATA_PIXEL_PERCENTAGE,
-        )
-        if date is None:
-            continue
-        else:
-            match_count += 1
-
-        if subset:
-            logger.info(f"Subsetting to {subset}")
-            pixels = [
-                part[:, subset[1] : subset[3], subset[0] : subset[2]] for part in pixels
-            ]
-
-        pixels = [part.compute() for part in pixels]
-
-        tiler(pixels, date, mgrs, bucket, localpath)
-
-        if match_count == DATES_PER_LOCATION:
-            break
-
-    if not match_count:
-        logger.info("No matching data found")
-
-
-if __name__ == "__main__":
-    main()
diff --git a/scripts/pipeline/landcover.py b/scripts/pipeline/landcover.py
deleted file mode 100644
index 68e4754a..00000000
--- a/scripts/pipeline/landcover.py
+++ /dev/null
@@ -1,201 +0,0 @@
-from pathlib import Path
-
-import click
-import fiona
-import geopandas
-import numpy
-import pandas
-import rasterio
-from fiona.crs import CRS
-from rasterio.windows import from_bounds
-from shapely.geometry import shape
-
-WGS84 = CRS.from_epsg(4326)
-NODATA = 0
-WATER = 80
-RANDOM_SEED = 23
-CLASSES = {
-    10: "Tree cover",
-    20: "Shrubland",
-    30: "Grassland",
-    40: "Cropland",
-    50: "Built-up",
-    60: "Bare / sparse vegetation",
-    70: "Snow and Ice",
-    80: "Permanent water bodies",
-    90: "Herbaceous wetland",
-    95: "Mangroves",
-    100: "Moss and lichen",
-}
-SCHEMA = {
-    "geometry": "MultiPolygon",
-    "properties": {
-        "name": "str",
-        "count": "int",
-        "Tree cover": "int",
-        "Shrubland": "int",
-        "Grassland": "int",
-        "Cropland": "int",
-        "Built-up": "int",
-        "Bare / sparse vegetation": "int",
-        "Snow and Ice": "int",
-        "Permanent water bodies": "int",
-        "Herbaceous wetland": "int",
-        "Mangroves": "int",
-        "Moss and lichen": "int",
-    },
-}
-
-
-@click.command()
-@click.option(
-    "--wd",
-    required=True,
-    type=str,
-)
-@click.option(
-    "--worldcover",
-    required=True,
-    type=str,
-)
-@click.option(
-    "--mgrs",
-    required=True,
-    type=str,
-)
-@click.option(
-    "--skip-stats",
-    is_flag=True,
-    required=False,
-    default=False,
-)
-def process(wd, worldcover, mgrs, skip_stats):
-    """
-    Run statistics and sampling.
-    """
-    if not skip_stats:
-        compute_stats(wd, worldcover, mgrs)
-    sample(wd)
-
-
-def compute_stats(wd, worldcover, mgrs):
-    """
-    Compute statistics of Worldcover data over MGRS tiles.
-    """
-    result = []
-    with rasterio.open(worldcover) as cover:
-        with fiona.open(mgrs, "r") as tiles:
-            assert cover.crs.to_epsg() == tiles.crs.to_epsg()
-            for tile in tiles:
-                print(tile.properties.get("Name"))
-                # Split polygons in parts, mgrs tiles at the dateline
-                # are split into two parts in the mgrs source file.
-                parts = shape(tile["geometry"]).geoms
-                pixels = []
-                for polygon in parts:
-                    bounds = from_bounds(*polygon.bounds, cover.transform)
-                    pixels.append(
-                        cover.read(
-                            1,
-                            window=bounds,
-                        ).ravel()
-                    )
-                pixels = numpy.hstack(pixels)
-
-                pixels = pixels[pixels != NODATA]
-
-                if not pixels.size:
-                    continue
-                elif numpy.all(pixels == WATER):
-                    continue
-                else:
-                    props = {}
-                    for key, classname in CLASSES.items():
-                        props[str(classname)] = int(numpy.sum(pixels == key))
-                    props["name"] = tile.properties.get("Name")
-                    props["count"] = int(len(numpy.unique(pixels)))
-
-                    result.append(
-                        {
-                            "geometry": dict(tile["geometry"]),
-                            "properties": props,
-                        }
-                    )
-
-    with fiona.open(
-        Path(wd, "mgrs_stats.fgb"),
-        "w",
-        driver="FlatGeobuf",
-        crs=WGS84,
-        schema=SCHEMA,
-    ) as colxn:
-        colxn.writerecords(result)
-
-
-def split_highest(data, column, size, pool=1000, seed=RANDOM_SEED):
-    """
-    Split highest values of a column from a dataframe.
-    """
-    data.sort_values(column, ascending=False, inplace=True)
-    return data[:pool].sample(size, random_state=seed)
-
-
-def percentages(data):
-    """
-    Normalize all numerical columns to percentages
-    """
-    data_num = data.select_dtypes(include="number")
-    data_norm = data_num.div(data_num.sum(axis=1), axis=0)
-    data[data_norm.columns] = data_norm
-
-    return data
-
-
-def sample(wd):
-    """
-    Sample the mgrs tiles based on landcover statistics.
-    """
-    data = geopandas.read_file(Path(wd, "mgrs_stats.fgb"))
-
-    data_norm = percentages(data.loc[:, data.columns != "count"])
-    data[data_norm.columns] = data_norm
-
-    diversity = split_highest(data, "count", 400, 2000)
-    urban = split_highest(data, "Built-up", 300, 300)
-    urban = split_highest(data, "Built-up", 1000, 1500)
-    wetland = split_highest(data, "Herbaceous wetland", 50, 500)
-    mangroves = split_highest(data, "Mangroves", 50, 500)
-    moss = split_highest(data, "Moss and lichen", 50, 500)
-    cropland = split_highest(data, "Cropland", 800, 3600)
-    trees = split_highest(data, "Tree cover", 150, 750)
-    shrubland = split_highest(data, "Shrubland", 100, 500)
-    grassland = split_highest(data, "Grassland", 200, 500)
-    bare = split_highest(data, "Bare / sparse vegetation", 50, 500)
-    snow = split_highest(data, "Snow and Ice", 25, 250)
-    water = split_highest(data, "Permanent water bodies", 50, 1000)
-
-    result = pandas.concat(
-        [
-            diversity,
-            urban,
-            wetland,
-            mangroves,
-            moss,
-            cropland,
-            trees,
-            shrubland,
-            grassland,
-            bare,
-            snow,
-            water,
-        ]
-    )
-
-    result = result.drop_duplicates(subset=["name"])
-    print(f"Found {len(result)} MGRS tiles")
-
-    result.to_file(Path(wd, "mgrs_sample.fgb", driver="FlatGeobuf"))
-
-
-if __name__ == "__main__":
-    process()
diff --git a/scripts/pipeline/landcover.sh b/scripts/pipeline/landcover.sh
deleted file mode 100644
index 8483a2e1..00000000
--- a/scripts/pipeline/landcover.sh
+++ /dev/null
@@ -1,27 +0,0 @@
-wd=/datadisk
-
-# Download Worldcover layer
-aws s3 sync s3://esa-worldcover/v200/2021/map $wd/esa-worldcover-v200-2021-map --no-sign-request
-
-# Download MGRS grid kml and convert to fgb
-curl -o $wd/S2A_OPER_GIP_TILPAR_MPC__20151209T095117_V20150622T000000_21000101T000000_B00.kml https://hls.gsfc.nasa.gov/wp-content/uploads/2016/03/S2A_OPER_GIP_TILPAR_MPC__20151209T095117_V20150622T000000_21000101T000000_B00.kml
-ogr2ogr \
-    -overwrite\
-    $wd/mgrs.fgb\
-    $wd/S2A_OPER_GIP_TILPAR_MPC__20151209T095117_V20150622T000000_21000101T000000_B00.kml \
-    -nlt multipolygon\
-    -sql "select Name, ExtractMultiPolygon(geometry) as geometry from Features"\
-    -dialect sqlite
-
-# Reduce resolution and merge
-counter=0
-for file in $wd/esa-worldcover-v200-2021-map/*.tif; do
-    counter=$((counter+1))
-    echo "$counter $file"
-    gdal_translate -ovr 3 $file $wd/tmp/tmp_$counter.tif
-done
-
-gdal_merge.py -o $wd/worldcover.tif $wd/tmp/*.tif
-
-# Intersect worldcover with mgrs and sample based on landcover
-python landcover.py --wd=$wd --worldcover=$wd/worldcover.tif --mgrs=$wd/mgrs.fgb
diff --git a/scripts/pipeline/tile.py b/scripts/pipeline/tile.py
deleted file mode 100644
index 6c3cf2a8..00000000
--- a/scripts/pipeline/tile.py
+++ /dev/null
@@ -1,146 +0,0 @@
-"""
-Tiling module for multi-dimensional datacubes.
-
-This script defines a tiling function for processing multi-dimensional imagery
-stacks into smaller tiles, while filtering out tiles with high cloud coverage
-or no-data pixels.
-
-It includes functions to filter tiles based on cloud coverage and no-data pixels,
-and a tiling function that generates smaller tiles from the input stack.
-"""
-
-import logging
-import subprocess
-import tempfile
-
-import numpy as np
-import rasterio
-import rioxarray  # noqa: F401
-import xarray as xr
-from rasterio.enums import ColorInterp
-
-NODATA = 0
-TILE_SIZE = 256
-PIXELS_PER_TILE = TILE_SIZE * TILE_SIZE
-BAD_PIXEL_MAX_PERCENTAGE = 0.3
-SCL_FILTER = [0, 1, 3, 8, 9, 10]
-VERSION = "04"
-
-logger = logging.getLogger("datacube")
-
-
-def filter_clouds_nodata(tile):
-    """
-    Filter tiles based on cloud coverage and no-data pixels.
-
-    Args:
-    - tile (xarray.Dataset): A subset of data representing a tile.
-
-    Returns:
-    - bool: True if the tile is approved, False if rejected.
-    """
-    # Check for nodata pixels
-    if int(tile.sel(band="B02").isin([NODATA]).sum()):
-        logger.debug("Too much no-data in B02")
-        return False
-
-    bands_to_check = ["vv", "vh", "dem"]
-    for band in bands_to_check:
-        if int(np.isnan(tile.sel(band=band)).sum()):
-            logger.debug(f"Too much no-data in {band}")
-            return False
-
-    # Check for cloud coverage
-    cloudy_pixel_count = int(tile.sel(band="SCL").isin(SCL_FILTER).sum())
-    if cloudy_pixel_count / PIXELS_PER_TILE >= BAD_PIXEL_MAX_PERCENTAGE:
-        logger.debug("Too much cloud coverage")
-        return False
-
-    return True  # If both conditions pass
-
-
-def tile_to_dir(stack, date, mgrs, bucket, dir):
-    """
-    Function to tile a multi-dimensional imagery stack while filtering out
-    tiles with high cloud coverage or no-data pixels.
-
-    Args:
-    - stack (xarray.Dataset): The input multi-dimensional imagery stack.
-    - date (str): Date string yyyy-mm-dd
-    - mgrs (str): MGRS Tile id
-    - bucket(str): AWS S3 bucket to write tiles to
-    """
-    logger.debug(f"Writing tempfiles to {dir}")
-
-    # Calculate the number of full tiles in x and y directions
-    num_x_tiles = stack[0].x.size // TILE_SIZE
-    num_y_tiles = stack[0].y.size // TILE_SIZE
-
-    counter = 0
-    # Iterate through each chunk of x and y dimensions and create tiles
-    for y_idx in range(num_y_tiles):
-        for x_idx in range(num_x_tiles):
-            # Calculate the start and end indices for x and y dimensions
-            # for the current tile
-            x_start = x_idx * TILE_SIZE
-            y_start = y_idx * TILE_SIZE
-            x_end = x_start + TILE_SIZE
-            y_end = y_start + TILE_SIZE
-
-            # Select the subset of data for the current tile
-            parts = [part[:, y_start:y_end, x_start:x_end] for part in stack]
-
-            # Only concat here to save memory, it converts S2 data to float
-            tile = xr.concat(parts, dim="band").rename("tile")
-
-            counter += 1
-            if counter % 250 == 0:
-                logger.info(f"Counted {counter} tiles")
-
-            if not filter_clouds_nodata(tile):
-                continue
-
-            tile = tile.drop_sel(band="SCL")
-
-            # Track band names and color interpretation
-            tile.attrs["long_name"] = [str(x.values) for x in tile.band]
-            color = [ColorInterp.blue, ColorInterp.green, ColorInterp.red] + [
-                ColorInterp.gray
-            ] * (len(tile.band) - 3)
-
-            # Write tile to tempdir
-            name = "{dir}/claytile_{mgrs}_{date}_v{version}_{counter}.tif".format(
-                dir=dir,
-                mgrs=mgrs,
-                date=date.replace("-", ""),
-                version=VERSION,
-                counter=str(counter).zfill(4),
-            )
-            tile.rio.to_raster(name, compress="deflate")
-
-            with rasterio.open(name, "r+") as rst:
-                rst.colorinterp = color
-                rst.update_tags(date=date)
-    if bucket:
-        logger.info(f"Syncing {dir} with s3://{bucket}/{VERSION}/{mgrs}/{date}")
-        subprocess.run(
-            [
-                "aws",
-                "s3",
-                "sync",
-                dir,
-                f"s3://{bucket}/{VERSION}/{mgrs}/{date}",
-                "--quiet",
-            ],
-            check=True,
-        )
-    else:
-        logger.info("No bucket specified, skipping S3 sync.")
-
-
-def tiler(stack, date, mgrs, bucket, dir):
-    if dir:
-        tile_to_dir(stack, date, mgrs, bucket, dir)
-    else:
-        with tempfile.TemporaryDirectory() as tmpdir:
-            tile_to_dir(stack, date, mgrs, bucket, tmpdir)
diff --git a/scripts/worldcover/embeddings_db.py b/scripts/worldcover/embeddings_db.py
deleted file mode 100644
index 1f106448..00000000
--- a/scripts/worldcover/embeddings_db.py
+++ /dev/null
@@ -1,58 +0,0 @@
-from pathlib import Path
-
-import geopandas as gpd
-import lancedb
-import matplotlib.pyplot as plt
-from skimage import io
-
-# Set working directory
-wd = "/home/usr/Desktop/"
-
-# To download the existing embeddings run aws s3 sync
-# aws s3 sync s3://clay-worldcover-embeddings /my/dir/clay-worldcover-embeddings
-
-vector_dir = Path(wd + "clay-worldcover-embeddings/v002/2021/")
-
-# Create new DB structure or open existing
-db = lancedb.connect(wd + "worldcoverembeddings_db")
-
-# Read all vector embeddings into a list
-data = []
-for strip in vector_dir.glob("*.gpq"):
-    print(strip)
-    tile_df = gpd.read_parquet(strip).to_crs("epsg:3857")
-
-    for _, row in tile_df.iterrows():
-        data.append(
-            {"vector": row["embeddings"], "year": 2021, "bbox": row.geometry.bounds}
-        )
-
-# Show table names
-db.table_names()
-
-# Drop existing table if exists
-db.drop_table("worldcover-2021-v001")
-
-# Create embeddings table and insert the vector data
-tbl = db.create_table("worldcover-2021-v001", data=data, mode="overwrite")
-
-
-# Visualize some image chips
-def plot(df, cols=10):
-    fig, axs = plt.subplots(1, cols, figsize=(20, 10))
-
-    for ax, (i, row) in zip(axs.flatten(), df.iterrows()):
-        bbox = row["bbox"]
-        url = f"https://services.terrascope.be/wms/v2?SERVICE=WMS&version=1.1.1&REQUEST=GetMap&layers=WORLDCOVER_2021_S2_TCC&BBOX={','.join([str(dat) for dat in bbox])}&SRS=EPSG:3857&FORMAT=image/png&WIDTH=512&HEIGHT=512"  # noqa: E501
-        image = io.imread(url)
-        ax.imshow(image)
-        ax.set_axis_off()
-
-    plt.tight_layout()
-    plt.show()
-
-
-# Select a vector by index, and search 10 similar pairs, and plot
-v = tbl.to_pandas()["vector"].values[10540]
-result = tbl.search(query=v).limit(5).to_pandas()
-plot(result, 5)
diff --git a/scripts/worldcover/run.py b/scripts/worldcover/run.py
deleted file mode 100755
index d1458e1b..00000000
--- a/scripts/worldcover/run.py
+++ /dev/null
@@ -1,251 +0,0 @@
-#!/usr/bin/env python3
-
-# import sys
-# sys.path.append("/home/tam/Documents/repos/model")
-
-import os
-import tempfile
-from math import floor
-
-import boto3
-import einops
-import geopandas as gpd
-import numpy
-import pyarrow as pa
-import rasterio
-import torch
-from rasterio.windows import Window
-from shapely import box
-from torchvision.transforms import v2
-
-from src.datamodule import ClayDataset
-from src.model_clay import CLAYModule
-
-YEAR = int(os.environ.get("YEAR", 2020))
-DATE = f"{YEAR}-06-01"
-TILE_SIZE = 12000
-CHIP_SIZE = 512
-E_W_INDEX_START = 67
-E_W_INDEX_END = 125
-N_S_INDEX_START = 24
-N_S_INDEX_END = 49
-YORIGIN = 50.0
-XORIGIN = -125.0
-PXSIZE = 8.333333333333333e-05
-
-RASTER_X_SIZE = (E_W_INDEX_END - E_W_INDEX_START) * TILE_SIZE
-RASTER_Y_SIZE = (N_S_INDEX_END - N_S_INDEX_START) * TILE_SIZE
-NODATA = 0
-CKPT_PATH = "s3://clay-model-ckpt/v0/mae_epoch-24_val-loss-0.46.ckpt"
-# CKPT_PATH = "https://huggingface.co/made-with-clay/Clay/resolve/main/Clay_v0.1_epoch-24_val-loss-0.46.ckpt"
-VERSION = "002"
-BUCKET = "clay-worldcover-embeddings"
-URL = "https://esa-worldcover-s2.s3.amazonaws.com/rgbnir/{year}/N{yidx}/ESA_WorldCover_10m_{year}_v{version}_N{yidx}W{xidx}_S2RGBNIR.tif"
-WC_VERSION_LOOKUP = {
-    2020: 100,
-    2021: 200,
-}
-
-
-MEAN = [
-    1369.03,  # red
-    1597.68,  # green
-    1741.10,  # blue
-    2858.43,  # nir
-]
-STD = [
-    2026.96,  # red
-    2011.88,  # green
-    2146.35,  # blue
-    2016.38,  # nir
-]
-
-grid = gpd.read_file(
-    # "/home/tam/Desktop/usa/esa_worldcover_grid_usa.fgb"
-    "https://clay-mgrs-samples.s3.amazonaws.com/esa_worldcover_grid_usa.fgb"
-)
-
-
-def tiles_and_windows(input: Window):
-    print("Input", input)
-    x_tile_index = E_W_INDEX_END - floor(input.col_off / TILE_SIZE)
-    x_local_off = input.col_off % TILE_SIZE
-    x_size = min(CHIP_SIZE, TILE_SIZE - x_local_off)
-    x_another = x_size < CHIP_SIZE
-
-    y_tile_index = N_S_INDEX_END - floor(input.row_off / TILE_SIZE)
-    y_local_off = input.row_off % TILE_SIZE
-    y_size = min(CHIP_SIZE, TILE_SIZE - y_local_off)
-    y_another = y_size < CHIP_SIZE
-
-    tile_id = f"N{y_tile_index}W{str(x_tile_index).zfill(3)}"
-    if tile_id not in grid.tile.values:
-        return
-
-    result = [
-        (
-            URL.format(
-                yidx=y_tile_index,
-                xidx=str(x_tile_index).zfill(3),
-                year=YEAR,
-                version=WC_VERSION_LOOKUP[YEAR],
-            ),
-            Window(x_local_off, y_local_off, x_size, y_size),
-        )
-    ]
-
-    if x_another:
-        result.append(
-            (
-                URL.format(
-                    yidx=y_tile_index,
-                    xidx=str(x_tile_index - 1).zfill(3),
-                    year=YEAR,
-                    version=WC_VERSION_LOOKUP[YEAR],
-                ),
-                Window(0, y_local_off, CHIP_SIZE - x_size, y_size),
-            )
-        )
-    if y_another:
-        result.append(
-            (
-                URL.format(
-                    yidx=y_tile_index - 1,
-                    xidx=str(x_tile_index).zfill(3),
-                    year=YEAR,
-                    version=WC_VERSION_LOOKUP[YEAR],
-                ),
-                Window(x_local_off, 0, x_size, CHIP_SIZE - y_size),
-            )
-        )
-    if x_another and y_another:
-        result.append(
-            (
-                URL.format(
-                    yidx=y_tile_index - 1,
-                    xidx=str(x_tile_index - 1).zfill(3),
-                    year=YEAR,
-                    version=WC_VERSION_LOOKUP[YEAR],
-                ),
-                Window(0, 0, CHIP_SIZE - x_size, CHIP_SIZE - y_size),
-            )
-        )
-
-    return result
-
-
-def make_batch(result):
-    pixels = []
-    for url, win in result:
-        with rasterio.open(url) as src:
-            data = src.read(window=win)
-            if NODATA in data:
-                return
-            pixels.append(data)
-            transform = src.window_transform(win)
-
-    if len(pixels) == 1:
-        pixels = pixels[0]
-    elif len(pixels) == 2:  # noqa: PLR2004
-        if pixels[0].shape[2] == CHIP_SIZE:
-            pixels = einops.pack(pixels, "b * w")[0]
-        else:
-            pixels = einops.pack(pixels, "b h *")[0]
-    else:
-        px1 = einops.pack(pixels[:2], "b w *")[0]
-        px2 = einops.pack(pixels[2:], "b w *")[0]
-        pixels = einops.pack((px1, px2), "b * w")[0]
-
-    assert pixels.shape == (4, CHIP_SIZE, CHIP_SIZE)
-
-    return {
-        "pixels": torch.as_tensor(data=[pixels], dtype=torch.float32).to(
-            rgb_model.device
-        ),
-        "latlon": torch.as_tensor(
-            data=[ds.normalize_latlon(transform[0], transform[3])]
-        ).to(rgb_model.device),
-        "timestep": torch.as_tensor(data=[ds.normalize_timestamp(f"{YEAR}-06-01")]).to(
-            rgb_model.device
-        ),
-    }
-
-
-index = int(os.environ.get("AWS_BATCH_JOB_ARRAY_INDEX", 2))
-
-# Setup model components
-tfm = v2.Compose([v2.Normalize(mean=MEAN, std=STD)])
-ds = ClayDataset(chips_path=[], transform=tfm)
-
-rgb_model = CLAYModule.load_from_checkpoint(
-    CKPT_PATH,
-    mask_ratio=0.0,
-    band_groups={"rgb": (0, 1, 2), "nir": (3,)},
-    strict=False,  # ignore the extra parameters in the checkpoint
-)
-
-xoff = index * CHIP_SIZE
-yoff = 0
-embeddings = []
-all_bounds = []
-while yoff < RASTER_Y_SIZE:
-    result = tiles_and_windows(Window(xoff, yoff, CHIP_SIZE, CHIP_SIZE))
-
-    if result is None:
-        yoff += CHIP_SIZE
-        continue
-
-    batch = make_batch(result)
-    if batch is None:
-        yoff += CHIP_SIZE
-        continue
-
-    (
-        unmasked_patches,
-        unmasked_indices,
-        masked_indices,
-        masked_matrix,
-    ) = rgb_model.model.encoder(batch)
-
-    embeddings.append(unmasked_patches.detach().cpu().numpy())
-    all_bounds.append(
-        (
-            XORIGIN + PXSIZE * xoff,
-            YORIGIN - PXSIZE * (yoff + CHIP_SIZE),
-            XORIGIN + PXSIZE * (xoff + CHIP_SIZE),
-            YORIGIN - PXSIZE * yoff,
-        )
-    )
-
-    yoff += CHIP_SIZE
-
-embeddings = numpy.vstack(embeddings)
-
-embeddings_mean = embeddings[:, :-2, :].mean(axis=1)
-
-print(f"Average embeddings have shape {embeddings_mean.shape}")
-
-gdf = gpd.GeoDataFrame(
-    data={
-        "embeddings": pa.FixedShapeTensorArray.from_numpy_ndarray(
-            numpy.ascontiguousarray(embeddings_mean)
-        ),
-    },
-    geometry=[box(*dat) for dat in all_bounds],  # This assumes same order
-    crs="EPSG:4326",
-)
-
-with tempfile.TemporaryDirectory() as tmp:
-    # tmp = "/home/tam/Desktop/wcctmp"
-
-    outpath = f"{tmp}/worldcover_embeddings_{YEAR}_{index}_v{VERSION}.gpq"
-    print(f"Uploading embeddings to {outpath}")
-
-    gdf.to_parquet(path=outpath, compression="ZSTD", schema_version="1.0.0")
-
-    s3_client = boto3.client("s3")
-    s3_client.upload_file(
-        outpath,
-        BUCKET,
-        f"v{VERSION}/{YEAR}/{os.path.basename(outpath)}",
-    )