Skip to content

darts_acquisition

darts_acquisition

Acquisition of data from various sources for the DARTS dataset.

__version__ module-attribute 

__version__ = importlib.metadata.version('darts-nextgen')

download_admin_files 

download_admin_files(admin_dir: pathlib.Path)

Download the admin files for the regions.

Files will be stored under [admin_dir]/adm1.shp and [admin_dir]/adm2.shp.

Parameters:

  • admin_dir (pathlib.Path) –

    The path to the admin files.

Source code in darts-acquisition/src/darts_acquisition/admin.py 
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
@stopwatch.f("Downloading admin files", printer=logger.debug)
def download_admin_files(admin_dir: Path):
    """Download the admin files for the regions.

    Files will be stored under [admin_dir]/adm1.shp and [admin_dir]/adm2.shp.

    Args:
        admin_dir (Path): The path to the admin files.

    """
    # Download the admin files
    admin_1_url = "https://github.com/wmgeolab/geoBoundaries/raw/main/releaseData/CGAZ/geoBoundariesCGAZ_ADM1.zip"
    admin_2_url = "https://github.com/wmgeolab/geoBoundaries/raw/main/releaseData/CGAZ/geoBoundariesCGAZ_ADM2.zip"

    admin_dir.mkdir(exist_ok=True, parents=True)

    logger.debug(f"Downloading {admin_1_url} to {admin_dir.resolve()}")
    _download_zip(admin_1_url, admin_dir)

    logger.debug(f"Downloading {admin_2_url} to {admin_dir.resolve()}")
    _download_zip(admin_2_url, admin_dir)

load_arcticdem 

load_arcticdem(
    geobox: odc.geo.geobox.GeoBox,
    data_dir: pathlib.Path | str,
    resolution: darts_acquisition.arcticdem.RESOLUTIONS,
    buffer: int = 0,
    persist: bool = True,
) -> xarray.Dataset

Load the ArcticDEM for the given geobox, fetch new data from the STAC server if necessary.

Parameters:

  • geobox (odc.geo.geobox.GeoBox) –

    The geobox for which the tile should be loaded.

  • data_dir (pathlib.Path | str) –

    The directory where the ArcticDEM data is stored.

  • resolution (typing.Literal[2, 10, 32]) –

    The resolution of the ArcticDEM data in m.

  • buffer (int, default: 0 ) –

    The buffer around the projected (epsg:3413) geobox in pixels. Defaults to 0.

  • persist (bool, default: True ) –

    If the data should be persisted in memory. If not, this will return a Dask backed Dataset. Defaults to True.

Returns:

  • xarray.Dataset –

    xr.Dataset: The ArcticDEM tile, with a buffer applied. Note: The buffer is applied in the arcticdem dataset's CRS, hence the orientation might be different. Final dataset is NOT matched to the reference CRS and resolution.

Warning

Geobox must be in a meter based CRS.

Usage

Since the API of the load_arcticdem is based on GeoBox, one can load a specific ROI based on an existing Xarray DataArray:

import xarray as xr
import odc.geo.xr

from darts_aquisition import load_arcticdem

# Assume "optical" is an already loaded s2 based dataarray

arcticdem = load_arcticdem(
    optical.odc.geobox,
    "/path/to/arcticdem-parent-directory",
    resolution=2,
    buffer=ceil(self.tpi_outer_radius / 2 * sqrt(2))
)

# Now we can for example match the resolution and extent of the optical data:
arcticdem = arcticdem.odc.reproject(optical.odc.geobox, resampling="cubic")

The buffer parameter is used to extend the region of interest by a certain amount of pixels. This comes handy when calculating e.g. the Topographic Position Index (TPI), which requires a buffer around the region of interest to remove edge effects.

Raises:

  • ValueError –

    If the resolution is not supported.

Source code in darts-acquisition/src/darts_acquisition/arcticdem.py 
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
@stopwatch.f("Loading ArcticDEM", printer=logger.debug, print_kwargs=["data_dir", "resolution", "buffer", "persist"])
def load_arcticdem(
    geobox: GeoBox, data_dir: Path | str, resolution: RESOLUTIONS, buffer: int = 0, persist: bool = True
) -> xr.Dataset:
    """Load the ArcticDEM for the given geobox, fetch new data from the STAC server if necessary.

    Args:
        geobox (GeoBox): The geobox for which the tile should be loaded.
        data_dir (Path | str): The directory where the ArcticDEM data is stored.
        resolution (Literal[2, 10, 32]): The resolution of the ArcticDEM data in m.
        buffer (int, optional): The buffer around the projected (epsg:3413) geobox in pixels. Defaults to 0.
        persist (bool, optional): If the data should be persisted in memory.
            If not, this will return a Dask backed Dataset. Defaults to True.

    Returns:
        xr.Dataset: The ArcticDEM tile, with a buffer applied.
            Note: The buffer is applied in the arcticdem dataset's CRS, hence the orientation might be different.
            Final dataset is NOT matched to the reference CRS and resolution.

    Warning:
        Geobox must be in a meter based CRS.

    Usage:
        Since the API of the `load_arcticdem` is based on GeoBox, one can load a specific ROI based on an existing Xarray DataArray:

        ```python
        import xarray as xr
        import odc.geo.xr

        from darts_aquisition import load_arcticdem

        # Assume "optical" is an already loaded s2 based dataarray

        arcticdem = load_arcticdem(
            optical.odc.geobox,
            "/path/to/arcticdem-parent-directory",
            resolution=2,
            buffer=ceil(self.tpi_outer_radius / 2 * sqrt(2))
        )

        # Now we can for example match the resolution and extent of the optical data:
        arcticdem = arcticdem.odc.reproject(optical.odc.geobox, resampling="cubic")
        ```

        The `buffer` parameter is used to extend the region of interest by a certain amount of pixels.
        This comes handy when calculating e.g. the Topographic Position Index (TPI), which requires a buffer around the region of interest to remove edge effects.

    Raises:
        ValueError: If the resolution is not supported.

    """  # noqa: E501
    odc.stac.configure_rio(cloud_defaults=True, aws={"aws_unsigned": True})

    assert ".icechunk" == data_dir.suffix, f"Data directory {data_dir} must have an .icechunk suffix!"

    match resolution:
        case 2:
            assert "2m" in data_dir.stem and "32m" not in data_dir.stem, (
                f"Data directory {data_dir} must have a '2m' in the name!"
            )
            accessor = smart_geocubes.ArcticDEM2m(data_dir)
        case 10:
            assert "10m" in data_dir.stem, f"Data directory {data_dir} must have a '10m' in the name!"
            accessor = smart_geocubes.ArcticDEM10m(data_dir)
        case 32:
            assert "32m" in data_dir.stem, f"Data directory {data_dir} must have a '32m' in the name!"
            accessor = smart_geocubes.ArcticDEM32m(data_dir)
        case _:
            raise ValueError(f"Resolution {resolution} not supported, only 2m, 10m and 32m are supported")

    accessor.assert_created()

    arcticdem = accessor.load(geobox, buffer=buffer, persist=persist)

    # Change dtype of the datamask to uint8 for later reproject_match
    arcticdem["arcticdem_data_mask"] = arcticdem.datamask.astype("uint8")

    # Clip values to -100, 3000 range (see docs about bands)
    arcticdem["dem"] = arcticdem["dem"].clip(-100, 3000)

    # Change dtype of arcticdem to float32 to save memory
    arcticdem = arcticdem.astype("float32")

    return arcticdem

load_planet_masks 

load_planet_masks(
    fpath: str | pathlib.Path,
) -> xarray.Dataset

Load the valid and quality data masks from a Planet scene.

Parameters:

  • fpath (str | pathlib.Path) –

    The file path to the Planet scene from which to derive the masks.

Raises:

  • FileNotFoundError –

    If no matching UDM-2 TIFF file is found in the specified path.

Returns:

  • xarray.Dataset –

    xr.Dataset: A merged xarray Dataset containing two data masks: - 'valid_data_mask': A mask indicating valid (1) and no data (0). - 'quality_data_mask': A mask indicating high quality (1) and low quality (0).

Source code in darts-acquisition/src/darts_acquisition/planet.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
@stopwatch.f("Loading Planet masks", printer=logger.debug)
def load_planet_masks(fpath: str | Path) -> xr.Dataset:
    """Load the valid and quality data masks from a Planet scene.

    Args:
        fpath (str | Path): The file path to the Planet scene from which to derive the masks.

    Raises:
        FileNotFoundError: If no matching UDM-2 TIFF file is found in the specified path.

    Returns:
        xr.Dataset: A merged xarray Dataset containing two data masks:
            - 'valid_data_mask': A mask indicating valid (1) and no data (0).
            - 'quality_data_mask': A mask indicating high quality (1) and low quality (0).

    """
    # Convert to Path object if a string is provided
    fpath = fpath if isinstance(fpath, Path) else Path(fpath)

    logger.debug(f"Loading data masks from {fpath.resolve()}")

    # Get imagepath
    udm_path = next(fpath.glob("*_udm2.tif"), None)
    if not udm_path:
        udm_path = next(fpath.glob("*_udm2_clip.tif"), None)
    if not udm_path:
        raise FileNotFoundError(f"No matching UDM-2 TIFF files found in {fpath.resolve()} (.glob('*_udm2.tif'))")

    # See udm classes here: https://developers.planet.com/docs/data/udm-2/
    da_udm = xr.open_dataarray(udm_path).astype("uint8")
    invalids = da_udm.sel(band=8).fillna(0) != 0
    low_quality = da_udm.sel(band=[2, 3, 4, 5, 6]).max(axis=0) == 1
    high_quality = ~low_quality & ~invalids
    qa_ds = (
        xr.where(high_quality, 2, 0)
        .where(~low_quality, 1)
        .where(~invalids, 0)
        .astype("uint8")
        .to_dataset(name="quality_data_mask")
        .drop_vars("band")
    )
    qa_ds["planet_udm"] = da_udm

    qa_ds["quality_data_mask"].attrs = {
        "data_source": "planet",
        "long_name": "Quality data mask",
        "description": "0 = Invalid, 1 = Low Quality, 2 = High Quality",
    }
    qa_ds["planet_udm"].attrs = {
        "data_source": "planet",
        "long_name": "Planet UDM",
        "description": "Usable Data Mask",
    }

    return qa_ds

load_planet_scene 

load_planet_scene(
    fpath: str | pathlib.Path,
) -> xarray.Dataset

Load a PlanetScope satellite GeoTIFF file and return it as an xarray datset.

Parameters:

  • fpath (str | pathlib.Path) –

    The path to the directory containing the TIFF files or a specific path to the TIFF file.

Returns:

Raises:

Source code in darts-acquisition/src/darts_acquisition/planet.py 
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
@stopwatch.f("Loading Planet scene", printer=logger.debug)
def load_planet_scene(fpath: str | Path) -> xr.Dataset:
    """Load a PlanetScope satellite GeoTIFF file and return it as an xarray datset.

    Args:
        fpath (str | Path): The path to the directory containing the TIFF files or a specific path to the TIFF file.

    Returns:
        xr.Dataset: The loaded dataset

    Raises:
        FileNotFoundError: If no matching TIFF file is found in the specified path.

    """
    # Convert to Path object if a string is provided
    fpath = fpath if isinstance(fpath, Path) else Path(fpath)

    # Check if the directory contains a PSOrthoTile or PSScene
    planet_type = parse_planet_type(fpath)
    logger.debug(f"Loading Planet PS {planet_type.capitalize()} from {fpath.resolve()}")

    # Get imagepath
    ps_image = next(fpath.glob("*_SR.tif"), None)
    if not ps_image:
        ps_image = next(fpath.glob("*_SR_clip.tif"), None)
    if not ps_image:
        raise FileNotFoundError(f"No matching TIFF files found in {fpath.resolve()} (.glob('*_SR.tif'))")

    ps_meta = next(fpath.glob("*_metadata.json"), None)
    if not ps_meta:
        raise FileNotFoundError(
            f"No matching metadata JSON files found in {fpath.resolve()} (.glob('*_metadata.json'))"
        )
    metadata = json.load(ps_meta.open())

    # Define band names and corresponding indices
    planet_da = xr.open_dataarray(ps_image)

    # Divide by 10000 to get reflectance between 0 and 1
    planet_da = planet_da.astype("float32") / 10000.0

    # Create a dataset with the bands
    bands = ["blue", "green", "red", "nir"]
    ds_planet = planet_da.assign_coords({"band": bands}).to_dataset(dim="band")
    for var in bands:
        ds_planet[var].attrs["long_name"] = f"PLANET {var.capitalize()}"
        ds_planet[var].attrs["units"] = "Reflectance"

    for var in ds_planet.data_vars:
        ds_planet[var].attrs["data_source"] = "planet"
        ds_planet[var].attrs["planet_type"] = planet_type

    # Add sun and elevation from metadata
    ds_planet.attrs["azimuth"] = metadata.get("sun_azimuth", float("nan"))
    ds_planet.attrs["elevation"] = metadata.get("sun_elevation", float("nan"))

    if planet_type == "scene":
        ds_planet.attrs["tile_id"] = fpath.stem
        ds_planet.attrs["planet_scene_id"] = fpath.stem
    elif planet_type == "orthotile":
        ds_planet.attrs["tile_id"] = f"{fpath.parent.stem}-{fpath.stem}"
        ds_planet.attrs["planet_orthotile_id"] = fpath.parent.stem
        ds_planet.attrs["planet_scene_id"] = fpath.stem

    return ds_planet

load_tcvis 

load_tcvis(
    geobox: odc.geo.geobox.GeoBox,
    data_dir: pathlib.Path | str,
    buffer: int = 0,
    persist: bool = True,
) -> xarray.Dataset

Load the TCVIS for the given geobox, fetch new data from GEE if necessary.

Parameters:

  • geobox (odc.geo.geobox.GeoBox) –

    The geobox to load the data for.

  • data_dir (pathlib.Path | str) –

    The directory to store the downloaded data for faster access for consecutive calls.

  • buffer (int, default: 0 ) –

    The buffer around the geobox in pixels. Defaults to 0.

  • persist (bool, default: True ) –

    If the data should be persisted in memory. If not, this will return a Dask backed Dataset. Defaults to True.

Returns:

Usage

Since the API of the load_tcvis is based on GeoBox, one can load a specific ROI based on an existing Xarray DataArray:

import xarray as xr
import odc.geo.xr

from darts_aquisition import load_tcvis

# Assume "optical" is an already loaded s2 based dataarray

tcvis = load_tcvis(
    optical.odc.geobox,
    "/path/to/tcvis-parent-directory",
)

# Now we can for example match the resolution and extent of the optical data:
tcvis = tcvis.odc.reproject(optical.odc.geobox, resampling="cubic")
Source code in darts-acquisition/src/darts_acquisition/tcvis.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
@stopwatch.f("Loading TCVIS", printer=logger.debug, print_kwargs=["data_dir", "buffer", "persist"])
def load_tcvis(
    geobox: GeoBox,
    data_dir: Path | str,
    buffer: int = 0,
    persist: bool = True,
) -> xr.Dataset:
    """Load the TCVIS for the given geobox, fetch new data from GEE if necessary.

    Args:
        geobox (GeoBox): The geobox to load the data for.
        data_dir (Path | str): The directory to store the downloaded data for faster access for consecutive calls.
        buffer (int, optional): The buffer around the geobox in pixels. Defaults to 0.
        persist (bool, optional): If the data should be persisted in memory.
            If not, this will return a Dask backed Dataset. Defaults to True.

    Returns:
        xr.Dataset: The TCVIS dataset.

    Usage:
        Since the API of the `load_tcvis` is based on GeoBox, one can load a specific ROI based on an existing Xarray DataArray:

        ```python
        import xarray as xr
        import odc.geo.xr

        from darts_aquisition import load_tcvis

        # Assume "optical" is an already loaded s2 based dataarray

        tcvis = load_tcvis(
            optical.odc.geobox,
            "/path/to/tcvis-parent-directory",
        )

        # Now we can for example match the resolution and extent of the optical data:
        tcvis = tcvis.odc.reproject(optical.odc.geobox, resampling="cubic")
        ```

    """  # noqa: E501
    assert ".icechunk" == data_dir.suffix, f"Data directory {data_dir} must have an .icechunk suffix!"
    accessor = smart_geocubes.TCTrend(data_dir, create_icechunk_storage=False)

    # We want to assume that the datacube is already created to be save in a multi-process environment
    accessor.assert_created()

    tcvis = accessor.load(geobox, buffer=buffer, persist=persist)

    # Rename to follow our conventions
    tcvis = tcvis.rename_vars(
        {
            "TCB_slope": "tc_brightness",
            "TCG_slope": "tc_greenness",
            "TCW_slope": "tc_wetness",
        }
    )

    return tcvis

parse_planet_type 

parse_planet_type(
    fpath: pathlib.Path,
) -> typing.Literal["orthotile", "scene"]

Parse the type of Planet data from the directory path.

Parameters:

  • fpath (pathlib.Path) –

    The directory path to the Planet data.

Returns:

  • typing.Literal['orthotile', 'scene'] –

    Literal["orthotile", "scene"]: The type of Planet data.

Raises:

  • ValueError –

    If the Planet data type cannot be parsed from the file path.

Source code in darts-acquisition/src/darts_acquisition/planet.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def parse_planet_type(fpath: Path) -> Literal["orthotile", "scene"]:
    """Parse the type of Planet data from the directory path.

    Args:
        fpath (Path): The directory path to the Planet data.

    Returns:
        Literal["orthotile", "scene"]: The type of Planet data.

    Raises:
        ValueError: If the Planet data type cannot be parsed from the file path.

    """
    # Cases for Scenes:
    # - YYYYMMDD_HHMMSS_NN_XXXX
    # - YYYYMMDD_HHMMSS_XXXX

    # Cases for Orthotiles:
    # NNNNNNN/NNNNNNN_NNNNNNN_YYYY-MM-DD_XXXX
    # NNNNNNN_NNNNNNN_YYYY-MM-DD_XXXX

    assert fpath.is_dir(), "fpath must be the parent directory!"

    ps_name_parts = fpath.stem.split("_")

    if len(ps_name_parts) == 3:
        # Must be scene or invalid
        date, time, ident = ps_name_parts
        if _is_valid_date(date, "%Y%m%d") and _is_valid_date(time, "%H%M%S") and len(ident) == 4:
            return "scene"

    if len(ps_name_parts) == 4:
        # Assume scene
        date, time, n, ident = ps_name_parts
        if _is_valid_date(date, "%Y%m%d") and _is_valid_date(time, "%H%M%S") and n.isdigit() and len(ident) == 4:
            return "scene"
        # Is not scene, assume orthotile
        chunkid, tileid, date, ident = ps_name_parts
        if chunkid.isdigit() and tileid.isdigit() and _is_valid_date(date, "%Y-%m-%d") and len(ident) == 4:
            return "orthotile"

    raise ValueError(
        f"Could not parse Planet data type from {fpath}."
        f"Expected a format of YYYYMMDD_HHMMSS_NN_XXXX or YYYYMMDD_HHMMSS_XXXX for scene, "
        "or NNNNNNN/NNNNNNN_NNNNNNN_YYYY-MM-DD_XXXX or NNNNNNN_NNNNNNN_YYYY-MM-DD_XXXX for orthotile."
        f"Got {fpath.stem} instead."
        "Please ensure that the parent directory of the file is used, instead of the file itself."
    )