# Managing Large Geospatial Arrays with TileDB

**URL:** <https://forum.tiledb.com/t/managing-large-geospatial-arrays-with-tiledb/624>\
**Category:** Uncategorized\
**Created:** [September 28, 2023, 2:05pm UTC](https://forum.tiledb.com/t/managing-large-geospatial-arrays-with-tiledb/624 "2023-09-28T14:05:13Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Baptiste\_Morel-Lab](https://yyz1.discourse-cdn.com/flex035/user_avatar/forum.tiledb.com/baptiste_morel-lab/32/293_2.png) [@Baptiste\_Morel-Lab](https://forum.tiledb.com/u/Baptiste_Morel-Lab)\
**Post date:** [September 28, 2023, 2:05pm UTC](https://forum.tiledb.com/t/managing-large-geospatial-arrays-with-tiledb/624/1 "2023-09-28T14:05:13Z")

</div>

I have a specific use case that I’d like to address with TileDB. I want to create a TileDB array that aggregates several rasters. My end goal is to be able to query the data using filters on dimensions or attributes, such as coordinates or a range of values. This could potentially represent a very large array. I’d like to be able to add data to this array progressively, and it could potentially cover the entire Earth in the future. This brings up the question of choosing an architecture that would optimally manage both write and read performance for this array.

We have already tried handling this using a sparse array. However, writing a single raster (with dimensions of 36000x36000 pixels) all at once seems to saturate our RAM, which has a capacity of 32GB. We then tried to distribute the writing by chunks using Dask, but the write operation becomes considerably slow (taking about ~15 minutes per image). We believe the slow performance is mainly due to the conversion operations to a 1D array.

We also experimented with a dense array and Dask’s `to_tiledb()` function, which appears to significantly improve performance. However, this introduces other challenges, such as the need to pre-create a dense, empty array. Additionally, there’s the limitation where we can’t load a subset with Dask and challenges with maintaining uniform types in dense array dimensions.

I’d appreciate any insights or recommendations on how to best approach this.

Here is one example of our code with sparse array:

```auto
with rasterio.open("./image.tiff", 'r') as src:
    data = src.read(1)
    print(data.shape)
    data = da.from_array(data, chunks='100 MiB')
    chunk_sizes = data.chunks
    print("Chunk sizes:", chunk_sizes)

    print("create linspace")
    longitudes = da.linspace(src.bounds.top, src.bounds.bottom, src.height, chunks=chunk_sizes[0])
    latitudes = da.linspace(src.bounds.left, src.bounds.right, src.width, chunks=chunk_sizes[1])
    lats, longs = da.meshgrid(latitudes, longitudes, indexing='ij')

    with tiledb.open("./sample.tdb", 'w') as array:
        for i in range(len(chunk_sizes[0])):
            for j in range(len(chunk_sizes[1])):
                chunk_data = data.blocks[i, j].compute().ravel()
                chunk_lats = lats.blocks[i, j].compute().ravel()
                chunk_longs = longs.blocks[i, j].compute().ravel()
                time_chunk = np.full(chunk_data.shape, 2021, dtype=np.int32)

                array[time_chunk, chunk_lats, chunk_longs] = {'data': chunk_data}

```

---

<div class="post-metadata">

**Author:** ![Nick\_Kules](https://yyz1.discourse-cdn.com/flex035/user_avatar/forum.tiledb.com/nick_kules/32/235_2.png) [@Nick\_Kules](https://forum.tiledb.com/u/Nick_Kules)\
**Post date:** [September 28, 2023, 7:05pm UTC](https://forum.tiledb.com/t/managing-large-geospatial-arrays-with-tiledb/624/2 "2023-09-28T19:05:29Z")

</div>

Hi Baptiste,

Great use case for TileDB. There were a few questions I wanted to clarify on your data and use case:

- In regards to aggregating rasters are you referring to performing a mosaic \ stitching operation?
- Is there a specific reason for selecting the sparse data model for this use?
  - If this is to accommodate gaps in data coverage, the dense model will actually still work fine as long as pixel size is always constant. Gaps will just be represented as nodata and will not have a disk space cost.
  - If it is in regards to dimension support, our recently added dimension labels could be an avenue to resolve this.

- Could you provide information on the array schema?
- Is this any sort of publicly available data that you can provide us to reproduce on our end?

We’re happy to hear more about this use case!

---

<div class="post-metadata">

**Author:** ![Baptiste\_Morel-Lab](https://yyz1.discourse-cdn.com/flex035/user_avatar/forum.tiledb.com/baptiste_morel-lab/32/293_2.png) [@Baptiste\_Morel-Lab](https://forum.tiledb.com/u/Baptiste_Morel-Lab)\
**Post date:** [September 29, 2023, 1:42pm UTC](https://forum.tiledb.com/t/managing-large-geospatial-arrays-with-tiledb/624/3 "2023-09-29T13:42:07Z")

</div>

Thank you for delving into our specific use case. Here’s the responses of your questions:

1. **Stitching vs. Mosaicking** : Currently, our focus is on stitching. Each of our rasters covers a distinct region without any overlap. When combined, these rasters span the entirety of the world.

2. **Sparse Array Choice** : Our decision to go with a sparse array wasn’t cemented on strong justifications. We initially chose this approach because:

3. **Array Schema** :

4. **Sample Data** :

Don’t hesitate if you have any other question and thanks again for your help.

---

<div class="post-metadata">

**Author:** ![Nick\_Kules](https://yyz1.discourse-cdn.com/flex035/user_avatar/forum.tiledb.com/nick_kules/32/235_2.png) [@Nick\_Kules](https://forum.tiledb.com/u/Nick_Kules)\
**Post date:** [September 29, 2023, 10:56pm UTC](https://forum.tiledb.com/t/managing-large-geospatial-arrays-with-tiledb/624/4 "2023-09-29T22:56:56Z")

</div>

Thanks Baptiste for the information and the data links. We will have a look at this data and provide some best practice examples for ingesting into TileDB with optimal performance in mind,

Once we have that information together would you be interested in jumping on a meeting to go over the suggestions? I believe we may have spoken before in the past. If you’re open to it I will send an email with meeting openings.