# Optimal TileDB Structure and Querying Procedure for Spaceborne Lidar (GEDI) Data

**URL:** <https://forum.tiledb.com/t/optimal-tiledb-structure-and-querying-procedure-for-spaceborne-lidar-gedi-data/737>\
**Category:** Uncategorized\
**Created:** [October 31, 2024, 5:24pm UTC](https://forum.tiledb.com/t/optimal-tiledb-structure-and-querying-procedure-for-spaceborne-lidar-gedi-data/737 "2024-10-31T17:24:58Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Simon\_Besnard](https://yyz1.discourse-cdn.com/flex035/user_avatar/forum.tiledb.com/simon_besnard/32/352_2.png) [@Simon\_Besnard](https://forum.tiledb.com/u/Simon_Besnard)\
**Post date:** [October 31, 2024, 5:24pm UTC](https://forum.tiledb.com/t/optimal-tiledb-structure-and-querying-procedure-for-spaceborne-lidar-gedi-data/737/1 "2024-10-31T17:24:58Z")

</div>

I am working on an optimal TileDB architecture and querying routine for spaceborne Lidar data from the GEDI mission. GEDI data, organized by orbit, includes variables at different levels of granularity:

- **Shot-level Variables** : Variables like `pai` (Plant Area Index) are stored at the shot level with `latitude`, `longitude`, and `time` dimensions.
- **Profile-level Variables** : Other variables, such as `rh` (relative height) and `pavd_z` (plant area volume density by height), include an additional vertical `profile_point` dimension to capture profiles within each shot.

I would like feedback on whether the following structure and reading/writing methods align with TileDB best practices.

### Current TileDB Structure and Methodology

1. **Data Structure and Array Setup** :

- **Two Array Structure** :
  - **Scalar Array** : Holds shot-level data using sparse indexing across `latitude`, `longitude`, and `time`.
  - **Profile Array** : Holds profile-level data with an additional `profile_point` dimension, allowing storage and retrieval of vertically indexed measurements for each shot.

- **Example Variables** :
  - Scalar variables: `pai`, `shot_number`
  - Profile variables: `rh`, `pavd_z`

1. **Writing Data** :

- Data is written in a sparse format for both arrays to support efficient indexing and range queries.
- The `profile_point` dimension in the profile array enables efficient storage and retrieval of profile information without redundancy.

1. **Reading Data** :

- Data retrieval uses `multi_index` with filtering latitude, longitude, and time range.
- Queries for profile data include the `profile_point` dimension for subsetting within each shot.

### Code Example for Reproducibility

The following code demonstrates the TileDB schema, write, and read processes I’m currently using.

#### Step 1: Initialize the Arrays

```auto

import tiledb
import numpy as np
import pandas as pd

# Configuration
lat_range = (-90.0, 90.0)
lon_range = (-180.0, 180.0)
time_range = (pd.Timestamp("2023-01-01").timestamp() * 1e6, pd.Timestamp("2023-12-31").timestamp() * 1e6)
max_profile_length = 101 # Profile levels per shot

# Scalar Array Schema
scalar_schema = tiledb.ArraySchema(
    domain=tiledb.Domain(
        tiledb.Dim(name="latitude", domain=lat_range, tile=0.1, dtype="float64"),
        tiledb.Dim(name="longitude", domain=lon_range, tile=0.1, dtype="float64"),
        tiledb.Dim(name="time", domain=time_range, tile=int(1e6), dtype="int64"),
    ),
    attrs=[tiledb.Attr(name="pai", dtype="float64")],
    sparse=True
)

# Profile Array Schema
profile_schema = tiledb.ArraySchema(
    domain=tiledb.Domain(
        tiledb.Dim(name="latitude", domain=lat_range, tile=0.1, dtype="float64"),
        tiledb.Dim(name="longitude", domain=lon_range, tile=0.1, dtype="float64"),
        tiledb.Dim(name="time", domain=time_range, tile=int(1e6), dtype="int64"),
        tiledb.Dim(name="profile_point", domain=(0, max_profile_length - 1), tile=1, dtype="int32")
    ),
    attrs=[
        tiledb.Attr(name="rh", dtype="float64"),
        tiledb.Attr(name="pavd_z", dtype="float64"),
        tiledb.Attr(name="shot_number", dtype="int64")
    ],
    sparse=True
)

# Create arrays (local URIs for example)
tiledb.Array.create("scalar_array", scalar_schema)
tiledb.Array.create("profile_array", profile_schema)

```

#### Step 2: Write Sample Data to Arrays

```auto

# Sample scalar data
scalar_data = {
    "latitude": np.random.uniform(lat_range[0], lat_range[1], 100),
    "longitude": np.random.uniform(lon_range[0], lon_range[1], 100),
    "time": np.random.randint(time_range[0], time_range[1], 100),
    "pai": np.random.rand(100)
}

# Write scalar data
with tiledb.open("scalar_array", mode="w") as array:
    array[scalar_data["latitude"], scalar_data["longitude"], scalar_data["time"]] = {"pai": scalar_data["pai"]}

# Sample profile data
profile_data = {
    "latitude": np.repeat(np.random.uniform(lat_range[0], lat_range[1], 10), max_profile_length),
    "longitude": np.repeat(np.random.uniform(lon_range[0], lon_range[1], 10), max_profile_length),
    "time": np.repeat(np.random.randint(time_range[0], time_range[1], 10), max_profile_length),
    "profile_point": np.tile(np.arange(max_profile_length), 10),
    "rh": np.random.rand(10 * max_profile_length),
    "pavd_z": np.random.rand(10 * max_profile_length),
    "shot_number": np.repeat(np.random.randint(1e10, 1e12, 10), max_profile_length)
}

# Write profile data
with tiledb.open("profile_array", mode="w") as array:
    array[profile_data["latitude"], profile_data["longitude"], profile_data["time"], profile_data["profile_point"]] = {
        "rh": profile_data["rh"],
        "pavd_z": profile_data["pavd_z"],
        "shot_number": profile_data["shot_number"]
    }

```

#### Step 3: Read Data Using `multi_index` Queries

```auto

def read_scalar_data(lat_min, lat_max, lon_min, lon_max, start_time, end_time, variables=["pai"]):
    with tiledb.open("scalar_array", mode="r") as array:
        query = array.query(attrs=variables)
        data = query.multi_index[lat_min:lat_max, lon_min:lon_max, start_time:end_time]
    print("Filtered Scalar Data:", data)

def read_profile_data(lat_min, lat_max, lon_min, lon_max, start_time, end_time, variables=["rh", "pavd_z", "shot_number"]):
    with tiledb.open("profile_array", mode="r") as array:
        query = array.query(attrs=variables)
        data = query.multi_index[lat_min:lat_max, lon_min:lon_max, start_time:end_time, :]
    print("Filtered Profile Data:", data)

# Example range values for querying
lat_min, lat_max = 5.0, 25.0
lon_min, lon_max = 25.0, 50.0
start_time = int(pd.Timestamp("2023-01-01").timestamp() * 1e6)
end_time = int(pd.Timestamp("2023-12-31").timestamp() * 1e6)

read_scalar_data(lat_min, lat_max, lon_min, lon_max, start_time, end_time)
read_profile_data(lat_min, lat_max, lon_min, lon_max, start_time, end_time)

```

### Questions for Optimization:

1. **Array Structure** :

- Given the two-array approach (scalar and profile arrays), would consolidating both data types into a single array with a `profile_point` dimension improve storage and access efficiency?

1. **Reading Efficiency** :

- The current querying uses `multi_index` on latitude, longitude, and time ranges. Are there alternative indexing or storage strategies in TileDB for enhancing query performance on such spatial-temporal data?

1. **Compression and Memory Management** :

- Any recommendations on compression and memory optimization when storing large profile data with extensive `profile_point` dimensions?

**Additional Question on Query Optimization for Large Spatial Areas** : Given that some queries will span large spatial regions and long timeframes, are there recommended strategies in TileDB to optimize read performance and memory usage for these cases? For instance:

- Are there specific tiling configurations or indexing techniques that work best for extensive spatial and temporal ranges?
- Would breaking down large spatial queries into smaller, more manageable regions or using a TileDB-native approach (e.g., subarray tiling) yield better performance?
- Any tips for balancing query speed and memory consumption when handling large, sparse datasets like GEDI?

Thanks a lot for you help!