diff --git a/.github/styles/config/vocabularies/docs/accept.txt b/.github/styles/config/vocabularies/docs/accept.txt index b856ecd..bfb955f 100644 --- a/.github/styles/config/vocabularies/docs/accept.txt +++ b/.github/styles/config/vocabularies/docs/accept.txt @@ -64,3 +64,5 @@ bool boolean (?i)modis reprojected +(?i)geospatial +(?i)cartesian diff --git a/api-reference/python/tilebox.datasets/Collection.find.mdx b/api-reference/python/tilebox.datasets/Collection.find.mdx index 654c4d4..2bf1103 100644 --- a/api-reference/python/tilebox.datasets/Collection.find.mdx +++ b/api-reference/python/tilebox.datasets/Collection.find.mdx @@ -42,6 +42,8 @@ data = collection.find( # check if a datapoint exists +from tilebox.datasets.sync.dataset import NotFoundError + try: collection.find( "0186d6b6-66cc-fcfd-91df-bbbff72499c3", diff --git a/assets/datasets/geometries/antimeridian-dark.png b/assets/datasets/geometries/antimeridian-dark.png new file mode 100644 index 0000000..1edc549 Binary files /dev/null and b/assets/datasets/geometries/antimeridian-dark.png differ diff --git a/assets/datasets/geometries/antimeridian-light.png b/assets/datasets/geometries/antimeridian-light.png new file mode 100644 index 0000000..706a8c3 Binary files /dev/null and b/assets/datasets/geometries/antimeridian-light.png differ diff --git a/assets/datasets/geometries/buggy-antimeridian-dark.png b/assets/datasets/geometries/buggy-antimeridian-dark.png new file mode 100644 index 0000000..0da95b2 Binary files /dev/null and b/assets/datasets/geometries/buggy-antimeridian-dark.png differ diff --git a/assets/datasets/geometries/buggy-antimeridian-light.png b/assets/datasets/geometries/buggy-antimeridian-light.png new file mode 100644 index 0000000..c615b93 Binary files /dev/null and b/assets/datasets/geometries/buggy-antimeridian-light.png differ diff --git a/assets/datasets/geometries/pole-covering-dark.png b/assets/datasets/geometries/pole-covering-dark.png new file mode 100644 index 0000000..2a48ba5 Binary files /dev/null and b/assets/datasets/geometries/pole-covering-dark.png differ diff --git a/assets/datasets/geometries/pole-covering-light.png b/assets/datasets/geometries/pole-covering-light.png new file mode 100644 index 0000000..330bd40 Binary files /dev/null and b/assets/datasets/geometries/pole-covering-light.png differ diff --git a/assets/datasets/geometries/polygon-both-poles-north-dark.png b/assets/datasets/geometries/polygon-both-poles-north-dark.png new file mode 100644 index 0000000..ec33678 Binary files /dev/null and b/assets/datasets/geometries/polygon-both-poles-north-dark.png differ diff --git a/assets/datasets/geometries/polygon-both-poles-north-light.png b/assets/datasets/geometries/polygon-both-poles-north-light.png new file mode 100644 index 0000000..8a9ce12 Binary files /dev/null and b/assets/datasets/geometries/polygon-both-poles-north-light.png differ diff --git a/assets/datasets/geometries/polygon-both-poles-south-dark.png b/assets/datasets/geometries/polygon-both-poles-south-dark.png new file mode 100644 index 0000000..19d269e Binary files /dev/null and b/assets/datasets/geometries/polygon-both-poles-south-dark.png differ diff --git a/assets/datasets/geometries/polygon-both-poles-south-light.png b/assets/datasets/geometries/polygon-both-poles-south-light.png new file mode 100644 index 0000000..3f12f5a Binary files /dev/null and b/assets/datasets/geometries/polygon-both-poles-south-light.png differ diff --git a/assets/datasets/geometries/winding-order-ccw.png b/assets/datasets/geometries/winding-order-ccw.png new file mode 100644 index 0000000..ce125e1 Binary files /dev/null and b/assets/datasets/geometries/winding-order-ccw.png differ diff --git a/assets/datasets/geometries/winding-order-cw.png b/assets/datasets/geometries/winding-order-cw.png new file mode 100644 index 0000000..d58d034 Binary files /dev/null and b/assets/datasets/geometries/winding-order-cw.png differ diff --git a/assets/datasets/queries/crs-cartesian-dark.png b/assets/datasets/queries/crs-cartesian-dark.png new file mode 100644 index 0000000..835f9a2 Binary files /dev/null and b/assets/datasets/queries/crs-cartesian-dark.png differ diff --git a/assets/datasets/queries/crs-cartesian-light.png b/assets/datasets/queries/crs-cartesian-light.png new file mode 100644 index 0000000..d3a81e2 Binary files /dev/null and b/assets/datasets/queries/crs-cartesian-light.png differ diff --git a/assets/datasets/queries/crs-spherical-dark.png b/assets/datasets/queries/crs-spherical-dark.png new file mode 100644 index 0000000..45fb72c Binary files /dev/null and b/assets/datasets/queries/crs-spherical-dark.png differ diff --git a/assets/datasets/queries/crs-spherical-light.png b/assets/datasets/queries/crs-spherical-light.png new file mode 100644 index 0000000..f100851 Binary files /dev/null and b/assets/datasets/queries/crs-spherical-light.png differ diff --git a/assets/datasets/queries/intersection-mode-contains-dark.png b/assets/datasets/queries/intersection-mode-contains-dark.png new file mode 100644 index 0000000..1564add Binary files /dev/null and b/assets/datasets/queries/intersection-mode-contains-dark.png differ diff --git a/assets/datasets/queries/intersection-mode-contains-light.png b/assets/datasets/queries/intersection-mode-contains-light.png new file mode 100644 index 0000000..8daf4ed Binary files /dev/null and b/assets/datasets/queries/intersection-mode-contains-light.png differ diff --git a/assets/datasets/queries/intersection-mode-intersects-dark.png b/assets/datasets/queries/intersection-mode-intersects-dark.png new file mode 100644 index 0000000..7e6bd84 Binary files /dev/null and b/assets/datasets/queries/intersection-mode-intersects-dark.png differ diff --git a/assets/datasets/queries/intersection-mode-intersects-light.png b/assets/datasets/queries/intersection-mode-intersects-light.png new file mode 100644 index 0000000..f417188 Binary files /dev/null and b/assets/datasets/queries/intersection-mode-intersects-light.png differ diff --git a/assets/geometries/antimeridian_buggy.png b/assets/geometries/antimeridian_buggy.png deleted file mode 100644 index 6cd3e38..0000000 Binary files a/assets/geometries/antimeridian_buggy.png and /dev/null differ diff --git a/assets/geometries/antimeridian_fixed.png b/assets/geometries/antimeridian_fixed.png deleted file mode 100644 index 05730dd..0000000 Binary files a/assets/geometries/antimeridian_fixed.png and /dev/null differ diff --git a/assets/geometries/granules.png b/assets/geometries/granules.png deleted file mode 100644 index 2cfd390..0000000 Binary files a/assets/geometries/granules.png and /dev/null differ diff --git a/assets/geometries/single.png b/assets/geometries/single.png deleted file mode 100644 index 6d8718c..0000000 Binary files a/assets/geometries/single.png and /dev/null differ diff --git a/assets/geometries/unary_union.png b/assets/geometries/unary_union.png deleted file mode 100644 index 52e49b6..0000000 Binary files a/assets/geometries/unary_union.png and /dev/null differ diff --git a/assets/guides/ingest/dataset-schema-dark.png b/assets/guides/ingest/dataset-schema-dark.png index e3a9106..9230410 100644 Binary files a/assets/guides/ingest/dataset-schema-dark.png and b/assets/guides/ingest/dataset-schema-dark.png differ diff --git a/assets/guides/ingest/dataset-schema-light.png b/assets/guides/ingest/dataset-schema-light.png index c871b55..94cb1ee 100644 Binary files a/assets/guides/ingest/dataset-schema-light.png and b/assets/guides/ingest/dataset-schema-light.png differ diff --git a/datasets/geometries.mdx b/datasets/geometries.mdx new file mode 100644 index 0000000..437c8ff --- /dev/null +++ b/datasets/geometries.mdx @@ -0,0 +1,415 @@ +--- +title: Working with Geometries +sidebarTitle: Geometries +description: Best practices for working with geometries in Tilebox. +icon: draw-polygon +--- + +Geometries are a common cause of friction in geospatial data processing. This is especially true when working with geometries that cross the antimeridian or cover a pole. +Tilebox is designed to take away most of the friction involved in this, but it's still recommended to follow the best practices for handling geometries outlined below. + +In doing so, you can ensure that no geometry related issues will arise even when interfacing with other libraries and tools that may not properly +support non-linearities in geometries. + +## Best Practices + +To ensure expected behavior when working with geometries, it's recommended to follow these best practices: + +1. Use [language specific geometry libraries](#geometry-libraries-and-common-formats) for creating, parsing and serializing geometries. +2. Polygon exteriors should always be defined in counter-clockwise [winding order](#winding-order). +3. Polygon holes (interior rings) should always be defined in clockwise [winding order](#winding-order). +4. Longitude values should always be within the `[-180, 180]` range, and latitude values should always be within the `[-90, 90]` range. +5. Geometries that cross the antimeridian should be [cut along the antimeridian](#antimeridian-cutting) into two parts—one for the eastern hemisphere and one for the western hemisphere. +7. Geometries that cover a pole need to be [defined in a specific way](#pole-coverings) to ensure correct handling. +8. When downloading geometries from external sources, such as from STAC Catalogues always make sure to verify that those assumptions are met to avoid unexpected behavior down the road. + +## Geometry vs Geography + +Some systems expose two similar, but different data types related to geometries: `Geometry` and `Geography`. The main difference lies in how edge interpolation along the antimeridian is handled. +`Geometry` represents a geometry in an arbitrary 2D cartesian coordinate system, and does not perform any edge interpolation, while `Geography` typically does wrap around the antimeridian by performing edge interpolation, such as converting cartesian coordinates to a 3D spherical coordinate system. +For `Geography` types, the `x` coordinate is typically limited to the `[-180, 180]` range, and the `y` coordinate is limited to the `[-90, 90]` range, while for `Geometry` types no such limitations are imposed. + +Tilebox does not make such a distinction between `Geometry` and `Geography` types. Instead, it provides a [query option to specify a coordinate reference system](/datasets/query/filter-by-location#coordinate-reference-system) to control whether geometry intersections and containment checks are performed in a 3D spherical coordinate system (which correctly handles antimeridian crossings) or in a standard 2D cartesian `lat/lon` coordinate system. + +## Terminology + +Get familiar with some key concepts related to geometries. + + + A great resource to learn more about these concepts is this [blog post by Tom MacWright](https://macwright.com/2015/03/23/geojson-second-bite). + + +**Point** + +A `Point` is a specific location on the Earth's surface, represented by a `longitude` and `latitude` coordinate. + +**Ring** + +A `Ring` is a collection of points that form a closed loop. The order of the points within the ring matter, since that defines the [winding order](#winding-order) of the ring, which is either clockwise or counter-clockwise. +Every ring should be explicitly closed, meaning that the first and last point should be the same. + +**Polygon** + +A `Polygon` is a collection of rings. The first ring is called the `exterior ring`, and represents the boundary of the polygon. Any other rings are called `interior rings`, and represent holes in the polygon. + +**MultiPolygon** + +A `MultiPolygon` is a collection of polygons. + +## Geometry libraries and common formats + +Geometries can be expressed in many different ways, and many formats exist for representing geometries. To handle these, the [Tilebox Client SDKs](/sdks/introduction) delegate geometry handling to external, well-known geometry libraries. + +| Client SDK | Geometry library used by Tilebox | +| ---------- | -------------------------------- | +| Python | [Shapely](https://shapely.readthedocs.io/en/stable/) | +| Go | [Orb](https://github.com/paulmach/orb) | + +Here is an example of how to define a `Polygon` geometry, covering the area of the state of Colorado using Tilebox Client SDKs. + + +```python Python +from shapely import Polygon + +area = Polygon( + ((-109.05, 41.00), (-109.045, 37.0), (-102.05, 37.0), (-102.05, 41.00), (-109.05, 41.00)), +) +``` +```go Go +// import "github.com/paulmach/orb" + +area := orb.Polygon{ + {{-109.05, 41.00}, {-109.045, 37.0}, {-102.05, 37.0}, {-102.05, 41.00}, {-109.05, 41.00}}, +} +``` + + +Through these libraries, Tilebox also supports some of the most common formats for representing geometries. + +### GeoJSON + +[GeoJSON](https://datatracker.ietf.org/doc/html/rfc7946) is a geospatial data interchange format based on JavaScript Object Notation (JSON). + +```json colorado.geojson +{ + "type": "Polygon", + "coordinates": [ + [ + [-109.05, 41.0], + [-109.045, 37.0], + [-102.05, 37.0], + [-102.05, 41.0], + [-109.05, 41.0] + ] + ] +} +``` + +Reading such a GeoJSON geometry can be achieved as follows. + + +```python Python +from shapely import from_geojson + +with open("colorado.geojson", "r") as f: + area = from_geojson(json.load(f)) +``` +```go Go +import ( + "github.com/paulmach/orb/geojson" + "os" +) + +func readGeoJSON() (orb.Geometry, error) { + geometryData, err := os.ReadFile("colorado.geojson") + if err != nil { + return nil, err + } + geojsonGeometry, err := geojson.UnmarshalGeometry(geometryData) + if err != nil { + return nil, err + } + return geojsonGeometry.Geometry(), nil +} +``` + + +### Well-Known Text (WKT) + +[Well-Known Text (WKT)](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry) is a text markup language for representing vector geometry objects. + + +```python Python +from shapely import from_wkt + +wkt = "POLYGON ((-109.05 41, -109.045 37, -102.05 37, -102.05 41, -109.05 41))" +area = from_wkt(wkt) +``` +```go Go +import ( + "github.com/paulmach/orb/encoding/wkt" +) + +func readWKT() (orb.Geometry, error) { + wktGeometry := "POLYGON ((-109.05 41, -109.045 37, -102.05 37, -102.05 41, -109.05 41))" + + return wkt.Unmarshal(wktGeometry) +} +``` + + +### Well-Known Binary (WKB) + +[Well-Known Binary (WKB)](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry#Well-known_binary) is a binary representation of vector geometry objects, a binary equivalent to WKT. + +Given a `colorado.wkb` file containing a binary encoding of the Colorado polygon, it can be read as follows. + + +```python Python +from shapely import from_wkb + +with open("colorado.wkb", "rb") as f: + area = from_wkb(f.read()) +``` +```go Go +import ( + "github.com/paulmach/orb/encoding/wkb" + "os" +) + +func readWKB() (orb.Geometry, error) { + binaryData, err := os.ReadFile("colorado.wkb") + if err != nil { + return nil, err + } + return wkb.Unmarshal(binaryData) +} +``` + + + + There is also an extended well known binary format (ewkb) that supports additional information such as specifying a spatial reference system (like EPSG:4326) in the encoded geometry. + Pythons `shapely` library supports that out of the box, and the `orb` library for Go supports it as well via the `github.com/paulmach/orb/encoding/ewkb` package. + + + +## Winding Order + +The winding order of a ring defines the orientation of the ring, and is either clockwise or counter-clockwise. The winding order of a ring is determined by the order of the points within the ring. + +Geometries in Tilebox follow the right hand rule defined by the [GeoJSON specification](https://datatracker.ietf.org/doc/html/rfc7946#section-3.1.6): + +{/* vale off */} +> A linear ring MUST follow the right-hand rule with respect to the area it bounds, i.e., exterior rings are counterclockwise, and holes are clockwise. +{/* vale on */} + +This means that the exterior ring of a polygon has a counter-clockwise winding order, while interior rings have a clockwise winding order. + +When this rule is not followed, querying may yield unexpected results, as can be seen in the example below. + +Take the following `Polygon` consisting of the same exterior ring but in different winding orders. + + +```plaintext Counter-clockwise +POLYGON ( + (-5 56, -6 29, 14 28, 23 55, -5 56) +) +``` +```plaintext Clockwise +POLYGON ( + (-5 56, 23 55, 14 28, -6 29, -5 56) +) +``` + + +This is how those two geometries would be interpreted on a sphere. + + + + Polygon over Europe with correct winding order + + + Polygon covering the whole globe with a hole over Europe due to incorrect winding order + + + + +Tilebox detects incorrect winding orders and automatically fixes them when ingesting data. + + +Since such unexpected behavior can cause issues and be hard to debug, Tilebox does detect incorrect winding orders and automatically fixes them when ingesting data. +This means that Geometries that really cover almost the whole globe except a small area **cannot** be ingested into Tilebox by simply specifying them as a `Polygon` with an exterior ring in clockwise winding order. + +Instead, such geometries should be defined as a `Polygon` consisting of two rings: +- an exterior ring covering the whole globe in counter-clockwise winding order +- an interior ring specifying the hole in clockwise winding order + +Such a definition, which Tilebox will interpret as intended, and also is valid according to the [GeoJSON specification](https://datatracker.ietf.org/doc/html/rfc7946#section-3.1.6), is shown below. + +```plaintext Polygon covering the whole globe with a hole over Europe +POLYGON ( + (-180 90, -180 -90, 180 -90, 180 90, -180 90), + (-5 56, 23 55, 14 28, -6 29, -5 56) +) +``` + +To verify the winding order of a geometry, you can use the following code snippets. + + +```python Python +from shapely import Polygon + +polygon = Polygon( + ((-109.05, 41.00), (-109.045, 37.0), (-102.05, 37.0), (-102.05, 41.00), (-109.05, 41.00)), +) +print(polygon.exterior.is_ccw) # True +``` +```go Go +import ( + "github.com/paulmach/orb" +) + +func isCounterClockwise() bool { + poly := orb.Polygon{ + {{-109.05, 41.00}, {-109.045, 37.0}, {-102.05, 37.0}, {-102.05, 41.00}, {-109.05, 41.00}}, + } + exterior := poly[0] + if exterior.Orientation() == orb.CCW { + return true + } + return false +} + +``` + + +## Antimeridian Crossings + +Geometries that cross the antimeridian are a common occurrence in satellite data, but can cause issues when not handled correctly. + +Take the following `Polygon` that crosses the antimeridian. + + + Polygon crossing the antimeridian + Polygon crossing the antimeridian + + +Below are a couple of different ways to express such a geometry. + +```plaintext Using exact longitude / latitude coordinates for each point +POLYGON ((173 5, -175 3, -172 20, 177 23, 173 5)) +``` + +While this is a valid representation, it causes issues for a lot of libraries that do calculations or visualization in a cartesian coordinate system, since the line from `lon=173` to `lon=-175` will be interpreted as a line that is `348` degrees long, crossing the null meridian and almost all the way around the globe. +Often times, visualizations of such geometries will look like in the image below. + + + Polygon crossing the antimeridian + Polygon crossing the antimeridian + + +To avoid such issues, some tools working in cartesian coordinate space sometimes extend the possible longitude range beyond to usual bounds of `[-180, 180]`. +Expressing the geometry in such a way, may look like this. + +```plaintext Extending the longitude range beyond 180 +POLYGON ((173 5, 185 3, 188 20, 177 23, 173 5)) +``` + +Or, the exact same area on the globe can also be expressed as a geometry by extending the longitude range below `-180`. + +```plaintext Extending the longitude range below -180 +POLYGON ((-187 5, -175 3, -172 20, -183 23, -187 5)) +``` + +While most visualization tools will probably handle such geometries correctly, special care needs to be taken when +working with such geometries when it comes to intersection and containment checks. + +The below code snippet illustrates this, by constructing a `Polygon` that covers the exact same area in the two different methods shown, and checking whether they intersect each other (which they should, since they represent the same area on the globe). + +```python Incorrect intersection check +from shapely import from_wkt + +a = from_wkt("POLYGON ((173 5, 185 3, 188 20, 177 23, 173 5))") +b = from_wkt("POLYGON ((-187 5, -175 3, -172 20, -183 23, -187 5))") + +print(a.intersects(b)) # False +``` + +### Antimeridian Cutting + +Additionally, none of the representations shown are valid according to the GeoJSON specification. The GeoJSON specification does offers a solution for this problem, though, which is [Antimeridian Cutting](https://datatracker.ietf.org/doc/html/rfc7946#section-3.1.9). +It suggests always cutting lines and polygons that cross the antimeridian into two parts—one for the eastern hemisphere and one for the western hemisphere. + +To achieve this, check out the [antimeridian](https://pypi.org/project/antimeridian/) python package, and an implementation of it in [Go](https://github.com/go-geospatial/antimeridian). + + + The `antimeridian` python package also is a great resource for [learning more about possible antimeridian issues](https://www.gadom.ski/antimeridian/latest/), how the [cutting algorithm](https://www.gadom.ski/antimeridian/latest/the-algorithm/) works and [edge cases](https://www.gadom.ski/antimeridian/latest/failure-modes/) to be aware of. + + +```plaintext Cutting the polygon along the antimeridian +MULTIPOLYGON ( + ((180 22.265994, 177 23, 173 5, 180 3.855573, 180 22.265994)), + ((-180 3.855573, -175 3, -172 20, -180 22.265994, -180 3.855573)) +) +``` + +By cutting geometries along the antimeridian: +- they conform to the GeoJSON specification +- there is no ambiguity about how to correctly interpret the geometry +- they are correctly handled by most visualization tools +- intersection and containment checks yield correct results no matter the [coordinate reference system](/datasets/query/filter-by-location#coordinate-reference-system) used +- all longitude values are within the valid `[-180, 180]` range, making it easier to work with them + + +## Pole Coverings + +Geometries that cover a pole can be especially tricky to handle correctly. No one algorithm that can handle all possible cases correctly, and different libraries and tools may handle them differently. +A lot of times, injecting prior knowledge, such as the fact that a geometry is supposed to cover a pole, can help to resolve such issues. For an example of this, check out the relevant section in +the [antimeridian documentation](https://www.gadom.ski/antimeridian/latest/failure-modes/#force-the-geometry-over-the-north-pole). + +Generally speaking though, there are two approaches that work well: + +### Approach 1: Cutting out a hole + +Define a geometry with an exterior ring that covers the whole globe. Then cut out a hole by defining an interior ring of the area that is not covered by the geometry. + + + This approach works especially well for geometries that cover both poles, since then the interior ring is guaranteed to not cover any of the poles. + + + + + Polygon covering both poles, viewing north pole covering + Polygon covering both poles, viewing north pole covering + + + Polygon covering both poles, viewing south pole covering + Polygon covering both poles, viewing south pole covering + + + +```plaintext Polygon covering both poles +POLYGON ( + (-180 90, -180 -90, 180 -90, 180 90, -180 90), + (-133 72, 153 77, 130 -64, -160 -66, -133 72) +) +``` + +### Approach 2: Splitting into multiple parts + +Another approach, that works well for circular caps covering a pole, involves cutting the geometry into multiple parts, but instead of splitting the geometry only at the antimeridian, also splitting it at the null meridian, the `90` and `-90` meridians. +For a circular cap covering the north pole, this results in four triangular parts, which can be combined into a `MultiPolygon`. Visualization libraries as well as intersection and containment checks done in cartesian coordinate space will then typically handle such a geometry correctly. + + + Polygon covering a pole + Polygon covering a pole + + +```plaintext Geometry of a circular cap covering the north pole +MULTIPOLYGON ( + ((-90 75, 0 75, 0 90, -90 90, -90 75)), + ((0 75, 90 75, 90 90, 0 90, 0 75)), + ((90 75, 180 75, 180 90, 90 90, 90 75)), + ((-180 75, -90 75, -90 90, -180 90, -180 75)) +) +``` diff --git a/datasets/ingest.mdx b/datasets/ingest.mdx index cf6be9a..3fd2b54 100644 --- a/datasets/ingest.mdx +++ b/datasets/ingest.mdx @@ -148,7 +148,7 @@ Measurements: [2025-03-28T11:44:23.000 UTC, 2025-03-28T11:45:19.000 UTC] (2 data #### xarray Dataset [`xarray.Dataset`](/sdks/python/xarray) is the default format in which Tilebox Datasets returns data when -[querying data](/datasets/query) from a collection. +[querying data](/datasets/query/querying-data) from a collection. Tilebox also supports it as input for ingestion. The example below shows how to construct an `xarray.Dataset` from scratch, that matches the schema of the `MyCustomDataset` dataset and can then be ingested into it. To learn more about `xarray.Dataset`, visit Tilebox dedicated [Xarray documentation page](/sdks/python/xarray). @@ -397,3 +397,9 @@ Through the usage of `xarray` and `pandas` you can also easily ingest existing d formats, such as CSV, [Parquet](https://parquet.apache.org/), [Feather](https://arrow.apache.org/docs/python/feather.html) and more. Check out the [Ingestion from common file formats](/guides/datasets/ingest-format) guide for examples of how to achieve this. + +## Geometries + +Ingesting Geometries can traditionally be a bit tricky, especially when working with geometries that cross the antimeridian or cover a pole. +Tilebox is designed to take away most of the friction involved in this, but it's still recommended to follow the [best practices for handling geometries](/datasets/geometries). + diff --git a/datasets/query.mdx b/datasets/query.mdx deleted file mode 100644 index b5a6205..0000000 --- a/datasets/query.mdx +++ /dev/null @@ -1,610 +0,0 @@ ---- -title: Querying data -sidebarTitle: Query -description: Learn how to query and load data from Tilebox datasets. -icon: server ---- - -Check out the examples below for common scenarios when loading data from collections. - - -```python Python -from tilebox.datasets import Client - -client = Client() -datasets = client.datasets() -collections = datasets.open_data.copernicus.sentinel1_sar.collections() -collection = collections["S1A_IW_RAW__0S"] -``` -```go Go -package main - -import ( - "context" - "log" - - "github.com/tilebox/tilebox-go/datasets/v1" -) - -func main() { - ctx := context.Background() - client := datasets.NewClient() - - dataset, err := client.Datasets.Get(ctx, "open_data.copernicus.sentinel1_sar") - if err != nil { - log.Fatalf("Failed to get dataset: %v", err) - } - - collection, err := client.Collections.Get(ctx, dataset.ID, "S1A_IW_RAW__0S") - if err != nil { - log.Fatalf("Failed to get collection: %v", err) - } -} -``` - - -To query data points from a dataset collection, use the [query](/api-reference/python/tilebox.datasets/Collection.query) method. It requires a `temporal_extent` parameter to specify the time or time interval for querying. - -## Filtering by time - -### Time interval - -To query data for a specific time interval, use a `tuple` in the form `(start, end)` as the `temporal_extent` parameter. Both `start` and `end` must be [TimeScalars](#time-scalars), which can be `datetime` objects or strings in ISO 8601 format. - - -```python Python -interval = ("2017-01-01", "2023-01-01") -data = collection.query(temporal_extent=interval, show_progress=True) -``` -```go Go -startDate := time.Date(2017, time.January, 1, 0, 0, 0, 0, time.UTC) -endDate := time.Date(2023, time.January, 1, 0, 0, 0, 0, time.UTC) -interval := query.NewTimeInterval(startDate, endDate) - -var datapoints []*v1.Sentinel1Sar -err = client.Datapoints.QueryInto(ctx, - []uuid.UUID{collection.ID}, - &datapoints, - datasets.WithTemporalExtent(interval), -) -if err != nil { - log.Fatalf("Failed to query datapoints: %v", err) -} - -log.Printf("Queried %d datapoints", len(datapoints)) -``` - - -Output - - -```plaintext Python - Size: 725MB -Dimensions: (time: 1109597, latlon: 2) -Coordinates: - ingestion_time (time) datetime64[ns] 9MB 2024-06-21T11:03:33.8524... - id (time) - - - The `show_progress` parameter is optional and can be used to display a [tqdm](https://tqdm.github.io/) progress bar while loading data. - - -A time interval specified as a tuple is interpreted as a half-closed interval. This means the start time is inclusive, and the end time is exclusive. -For instance, using an end time of `2023-01-01` includes data points up to `2022-12-31 23:59:59.999`, but excludes those from `2023-01-01 00:00:00.000`. -This behavior mimics the Python `range` function and is useful for chaining time intervals. - - -```python Python -import xarray as xr - -data = [] -for year in [2017, 2018, 2019, 2020, 2021, 2022]: - interval = (f"{year}-01-01", f"{year + 1}-01-01") - data.append(collection.query(temporal_extent=interval, show_progress=True)) - -# Concatenate the data into a single dataset, which is equivalent -# to the result of the single request in the code example above. -data = xr.concat(data, dim="time") -``` -```go Go -var datapoints []*v1.Sentinel1Sar - -for year := 2017; year <= 2022; year++ { - startDate := time.Date(year, time.January, 1, 0, 0, 0, 0, time.UTC) - interval := query.NewTimeInterval(startDate, startDate.AddDate(1, 0, 0)) - - var partialDatapoints []*v1.Sentinel1Sar - err = client.Datapoints.QueryInto(ctx, - []uuid.UUID{collection.ID}, - &partialDatapoints, - datasets.WithTemporalExtent(interval), - ) - if err != nil { - log.Fatalf("Failed to query datapoints: %v", err) - } - - // Concatenate the data into a single dataset, which is equivalent - // to the result of the single request in the code example above. - datapoints = append(datapoints, partialDatapoints...) -} -``` - - -Above example demonstrates how to split a large time interval into smaller chunks while loading data in separate requests. Typically, this is not necessary as the datasets client auto-paginates large intervals. - -### Endpoint inclusivity - -For greater control over inclusivity of start and end times, you can use the `TimeInterval` dataclass instead of a tuple of two [TimeScalars](#time-scalars). This class allows you to specify the `start` and `end` times, as well as their inclusivity. Here's an example of creating equivalent `TimeInterval` objects in two different ways. - - -```python Python -from datetime import datetime -from tilebox.datasets.data import TimeInterval - -interval1 = TimeInterval( - datetime(2017, 1, 1), datetime(2023, 1, 1), - end_inclusive=False -) -interval2 = TimeInterval( - # python datetime granularity is in milliseconds - datetime(2017, 1, 1), datetime(2022, 12, 31, 23, 59, 59, 999999), - end_inclusive=True -) - -print("Inclusivity is indicated by interval notation: ( and [") -print(interval1) -print(interval2) -print(f"They are equivalent: {interval1 == interval2}") -print(interval2.to_half_open()) - -# Query data for a time interval -data = collection.query(temporal_extent=interval1, show_progress=True) -``` -```go Go -interval1 := query.TimeInterval{ - Start: time.Date(2017, time.January, 1, 0, 0, 0, 0, time.UTC), - End: time.Date(2023, time.January, 1, 0, 0, 0, 0, time.UTC), - EndInclusive: false, -} - -interval2 := query.TimeInterval{ - Start: time.Date(2017, time.January, 1, 0, 0, 0, 0, time.UTC), - End: time.Date(2022, time.December, 31, 23, 59, 59, 999999999, time.UTC), - EndInclusive: true, -} - -log.Println("Inclusivity is indicated by interval notation: ( and [") -log.Println(interval1.String()) -log.Println(interval2.String()) -log.Println("They are equivalent:", interval1.Equal(&interval2)) -log.Println(interval2.ToHalfOpen().String()) - -// Query data for a time interval -var datapoints []*v1.Sentinel1Sar -err = client.Datapoints.QueryInto(ctx, - []uuid.UUID{collection.ID}, - &datapoints, - datasets.WithTemporalExtent(interval1), -) -``` - - -```plaintext Output -Inclusivity is indicated by interval notation: ( and [ -[2017-01-01T00:00:00.000 UTC, 2023-01-01T00:00:00.000 UTC) -[2017-01-01T00:00:00.000 UTC, 2022-12-31T23:59:59.999 UTC] -They are equivalent: True -[2017-01-01T00:00:00.000 UTC, 2023-01-01T00:00:00.000 UTC) -``` - -### Time scalars - -You can load all datapoints linked to a specific timestamp by specifying a `TimeScalar` as the time query argument. A `TimeScalar` can be a `datetime` object or a string in ISO 8601 format. When passed to the `load` method, it retrieves all data points matching exactly that specified time, with a millisecond precision. - - - A collection may contain multiple datapoints for one millisecond, so multiple data points could still be returned. If you want to fetch only a single data point, [query the collection by id](#loading-a-data-point-by-id) instead. - - -Here's how to query a data point at a specific millisecond from a [collection](/datasets/concepts/collections). - - - ```python Python - data = collection.query(temporal_extent="2024-08-01 00:00:01.362") - print(data) - ``` - ```go Go - temporalExtent := query.NewPointInTime(time.Date(2024, time.August, 1, 0, 0, 1, 362000000, time.UTC)) - - var datapoints []*v1.Sentinel1Sar - err = client.Datapoints.QueryInto(ctx, - []uuid.UUID{collection.ID}, &datapoints, - datasets.WithTemporalExtent(temporalExtent), - ) - - log.Printf("Queried %d datapoints", len(datapoints)) - log.Printf("First datapoint time: %s", datapoints[0].GetTime().AsTime()) - ``` - - -Output - - -```plaintext Python - Size: 721B -Dimensions: (time: 1, latlon: 2) -Coordinates: - ingestion_time (time) datetime64[ns] 8B 2024-08-01T08:53:08.450499 - id (time) - - - Tilebox uses millisecond precision for timestamps. To load all data points for a specific second, it's a [time interval](/datasets/query#time-interval) request. Refer to the examples below for details. - - -The output of the `query` method is an `xarray.Dataset` object. To learn more about Xarray, visit the dedicated [Xarray page](/sdks/python/xarray). - -### Time iterables (Python only) - -You can specify a time interval by using an iterable of `TimeScalar`s as the `temporal_extent` parameter. This is especially useful when you want to use the output of a previous `query` call as input for another query. Here's how that works. - - - ```python Python - interval = ("2017-01-01", "2023-01-01") - found_datapoints = collection.query(temporal_extent=interval, skip_data=True) - - first_50_data_points = collection.query(temporal_extent=found_datapoints.time[:50]) - print(first_50_data_points) - ``` - - -```plaintext Output - Size: 33kB -Dimensions: (time: 50, latlon: 2) -Coordinates: - ingestion_time (time) datetime64[ns] 400B 2024-06-21T11:03:33.852... - id (time) -```python Python -from datetime import datetime -import pytz - -# Tokyo has a UTC+9 hours offset, so this is the same as -# 2017-01-01 02:45:25.679 UTC -tokyo_time = pytz.timezone('Asia/Tokyo').localize( - datetime(2017, 1, 1, 11, 45, 25, 679000) -) -print(tokyo_time) -data = collection.query(temporal_extent=tokyo_time) -print(data) -``` -```go Go -// Tokyo has a UTC+9 hours offset, so this is the same as -// 2017-01-01 02:45:25.679 UTC -location, _ := time.LoadLocation("Asia/Tokyo") -tokyoTime := query.NewPointInTime(time.Date(2017, 1, 1, 11, 45, 25, 679000000, location)) -log.Println(tokyoTime) - -var datapoints []*v1.Sentinel1Sar -err = client.Datapoints.QueryInto(ctx, - []uuid.UUID{collection.ID}, &datapoints, - datasets.WithTemporalExtent(tokyoTime), -) -if err != nil { - log.Fatalf("Failed to query datapoints: %v", err) -} - -log.Printf("Queried %d datapoints", len(datapoints)) -// time is in UTC since API always returns UTC timestamps -log.Printf("First datapoint time: %s", datapoints[0].GetTime().AsTime()) -``` - - -Output - - -```plaintext Python -2017-01-01 11:45:25.679000+09:00 - Size: 725B -Dimensions: (time: 1, latlon: 2) -Coordinates: - ingestion_time (time) datetime64[ns] 8B 2024-06-21T11:03:33.852435 - id (time) - -## Filtering by area of interest - -[Spatio-temporal](/datasets/types/spatiotemporal) also come with spatial filtering capabilities. When querying, you can specify a time interval, and additionally also specify a bounding box or a polygon as an area of interest to filter by. - -Here is how to query Sentinel-2 `S2A_S2MSI2A` data over Colorado for April 2025. - - -```python Python -from shapely import MultiPolygon -from tilebox.datasets import Client - -area = MultiPolygon( - [ - (((-109.10, 40.98), (-102.01, 40.95), (-102.06, 36.82), (-109.06, 36.96), (-109.10, 40.98)),), - ] -) - -client = Client() -datasets = client.datasets() -sentinel2_msi = datasets.open_data.copernicus.sentinel2_msi -data = sentinel2_msi.collection("S2A_S2MSI2A").query( - temporal_extent=("2025-04-01", "2025-05-02"), - spatial_extent=area, - show_progress=True, -) -``` -```go Go -ctx := context.Background() -client := datasets.NewClient() - -dataset, err := client.Datasets.Get(ctx, "open_data.copernicus.sentinel2_msi") -if err != nil { - log.Fatalf("Failed to get dataset: %v", err) -} - -collection, err := client.Collections.Get(ctx, dataset.ID, "S2A_S2MSI2A") -if err != nil { - log.Fatalf("Failed to get collection: %v", err) -} - -startDate := time.Date(2025, 4, 1, 0, 0, 0, 0, time.UTC) -endDate := time.Date(2025, 5, 2, 0, 0, 0, 0, time.UTC) -timeInterval := query.NewTimeInterval(startDate, endDate) -area := orb.MultiPolygon{ - { - {{-109.10, 40.98}, {-102.01, 40.95}, {-102.06, 36.82}, {-109.06, 36.96}, {-109.10, 40.98}}, - }, -} - -var datapoints []*v1.Sentinel2Msi -err = client.Datapoints.QueryInto(ctx, - []uuid.UUID{collection.ID}, &datapoints, - datasets.WithTemporalExtent(timeInterval), - datasets.WithSpatialExtent(area), -) -if err != nil { - log.Fatalf("Failed to query datapoints: %v", err) -} -``` - - -## Skipping data fields - -Sometimes, only the ID or timestamp associated with a datapoint is required. In this case, loading the full data fields for each datapoint is not necessary and can be avoided by -setting the `skip_data` parameter to `True`. - -For example, when only checking how many datapoints exist in a given time interval, you can use `skip_data=True` to avoid loading the data fields. - - - ```python Python - interval = ("2023-01-01", "2023-02-01") - data = collection.query(temporal_extent=interval, skip_data=True) - print(f"Found {data.sizes['time']} data points.") - ``` -```go Go -startDate := time.Date(2023, time.January, 1, 0, 0, 0, 0, time.UTC) -endDate := time.Date(2023, time.February, 1, 0, 0, 0, 0, time.UTC) -interval := query.NewTimeInterval(startDate, endDate) - -var datapoints []*v1.Sentinel1Sar -err = client.Datapoints.QueryInto(ctx, - []uuid.UUID{collection.ID}, &datapoints, - datasets.WithTemporalExtent(temporalExtent), - datasets.WithSkipData(), -) -if err != nil { - log.Fatalf("Failed to query datapoints: %v", err) -} - -log.Printf("Queried %d datapoints", len(datapoints)) -log.Printf("First datapoint time: %s", datapoints[0].GetTime().AsTime()) -``` - - -Output - - -```plaintext Python - Size: 160B -Dimensions: (time: 1) -Coordinates: - ingestion_time (time) datetime64[ns] 8B 2024-08-01T08:53:08.450499 - id (time) - -## Empty response - -The `query` method always returns an `xarray.Dataset` object, even if there are no data points for the specified query. In such cases, the returned dataset will be empty, but no error will be raised. - - - ```python Python - time_with_no_data_points = "1997-02-06 10:21:00" - data = collection.query(temporal_extent=time_with_no_data_points) - print(data) - ``` -```go Go -timeWithNoDatapoints := query.NewPointInTime(time.Date(1997, time.February, 6, 10, 21, 0, 0, time.UTC)) - -var datapoints []*v1.Sentinel1Sar -err = client.Datapoints.QueryInto(ctx, - []uuid.UUID{collection.ID}, &datapoints, - datasets.WithTemporalExtent(timeWithNoDatapoints), -) -if err != nil { - log.Fatalf("Failed to query datapoints: %v", err) -} - -log.Printf("Queried %d datapoints", len(datapoints)) -``` - - -Output - - -```plaintext Python - Size: 0B -Dimensions: () -Data variables: - *empty* -``` -```plaintext Go -Queried 0 datapoints -``` - - -## By datapoint ID - -If you know the ID of the data point you want to load, you can use [find](/api-reference/python/tilebox.datasets/Collection.find). - -This method always returns a single data point or raises an exception if no data point with the specified ID exists. - - - ```python Python - datapoint_id = "01916d89-ba23-64c9-e383-3152644bcbde" - datapoint = collection.find(datapoint_id) - print(datapoint) - ``` -```go Go -datapointID := uuid.MustParse("01910b3c-8552-424d-e116-81d0c3402ccc") - -var datapoint v1.Sentinel1Sar -err = client.Datapoints.GetInto(ctx, - []uuid.UUID{collection.ID}, datapointID, &datapoint, -) -if err != nil { - log.Fatalf("Failed to query datapoint: %v", err) -} - -fmt.Println(protojson.Format(&datapoint)) -``` - - -Output - - -```plaintext Python - Size: 725B -Dimensions: (latlon: 2) -Coordinates: - ingestion_time datetime64[ns] 8B 2024-08-20T05:53:08.600528 - id - - - You can also set the `skip_data` parameter when calling `find` to query only the required fields of the data point, same as for `query`. - - -## Automatic pagination - -Querying large time intervals can return a large number of data points. -Tilebox automatically handles pagination for you by sending paginated requests to the server. - - -When using the python SDK in an interactive notebook environment, you can additionally also display a -progress bar to keep track of the progress of the query by setting the `show_progress` parameter to `True`. - diff --git a/datasets/query/filter-by-id.mdx b/datasets/query/filter-by-id.mdx new file mode 100644 index 0000000..8ed2517 --- /dev/null +++ b/datasets/query/filter-by-id.mdx @@ -0,0 +1,134 @@ +--- +title: Querying individual datapoints by ID +sidebarTitle: Filter by ID +description: Learn how to run query for a specific datapoint by its unique ID. +icon: fingerprint +--- + +If you already know the ID of the data point you want to query, you can query it directly, using +[Collection.find](/api-reference/python/tilebox.datasets/Collection.find) in Python or [Datapoints.GetInto](/api-reference/go/datasets/Datapoints.GetInto) in Go. + + + Check out [selecting a collection](/datasets/query/querying-data#selecting-a-collection) to learn how to get a collection object + to query from. + + + +```python Python +datapoint_id = "0197a491-1520-102f-48f4-f087d6ef8603" +datapoint = collection.find(datapoint_id) +print(datapoint) +``` +```go Go +datapointID := uuid.MustParse("0197a491-1520-102f-48f4-f087d6ef8603") + +var datapoint v1.Sentinel2Msi +err = client.Datapoints.GetInto(ctx, + []uuid.UUID{collection.ID}, datapointID, &datapoint, +) +if err != nil { + log.Fatalf("Failed to query datapoint: %v", err) +} + +fmt.Println(protojson.Format(&datapoint)) +``` + + +Output + + +```plaintext Python + Size: 443B +Dimensions: () +Coordinates: + time datetime64[ns] 8B 2025-06-25T00:51:01.024000 +Data variables: (12/23) + id + + + You can also set the `skip_data` parameter when calling `find` to query only the required fields of the data point, same as for `query`. + + +## Checking if a datapoint exists + +`find` returns an error if the specified datapoint does not exist. You can use this to check if a datapoint exists or not. + + +```python Python +from tilebox.datasets.sync.dataset import NotFoundError + +datapoint_id = "0197a47f-a830-1160-6df5-61ac723dae17" # doesn't exist + +try: + collection.find(datapoint_id) + exists = True +except NotFoundError: + exists = False +``` +```go Go +datapointID := uuid.MustParse("0197a47f-a830-1160-6df5-61ac723dae17") + +exists := True +var datapoint v1.Sentinel2Msi +err = client.Datapoints.GetInto(ctx, + []uuid.UUID{collection.ID}, datapointID, &datapoint, +) +if err != nil { + if connect.CodeOf(err) == connect.CodeNotFound { + exists = false + } else { + log.Fatalf("Failed to query datapoint: %v", err) + } +} +``` + diff --git a/datasets/query/filter-by-location.mdx b/datasets/query/filter-by-location.mdx new file mode 100644 index 0000000..ab9d3b4 --- /dev/null +++ b/datasets/query/filter-by-location.mdx @@ -0,0 +1,301 @@ +--- +title: Filtering by a location +sidebarTitle: Filter by location +description: Learn how to filter your query results by an area of interest specified as certain geographical extent. +icon: globe-pointer +--- + + + Check out the [best practices for handling geometries](/datasets/geometries) to learn more about the different aspects to consider when working with geometries, including antimeridian crossings and pole coverings. + + +When querying, you can specify arbitrary geometries as an area of interest to filter by. Tilebox currently supports +`Polygon` and `MultiPolygon` geometries as query filters. + +## Filtering by an area of interest + +To filter by an area of interest, use either a `Polygon` or `MultiPolygon` geometry as the spatial extent parameter. + +Here is how to query Sentinel-2 `S2A_S2MSI2A` data over Colorado for a certain day in April 2025. + + +```python Python +from shapely import Polygon +from tilebox.datasets import Client + +area = Polygon( # area roughly covering the state of Colorado + ((-109.05, 41.00), (-109.045, 37.0), (-102.05, 37.0), (-102.05, 41.00), (-109.05, 41.00)), +) + +client = Client() +sentinel2_msi = client.dataset("open_data.copernicus.sentinel2_msi") +collection = sentinel2_msi.collection("S2A_S2MSI2A") +data = collection.query( + temporal_extent=("2025-04-02", "2025-04-03"), + spatial_extent=area, +) +``` +```go Go +startDate := time.Date(2025, 4, 2, 0, 0, 0, 0, time.UTC) +endDate := time.Date(2025, 4, 3, 0, 0, 0, 0, time.UTC) +area := orb.Polygon{ // area roughly covering the state of Colorado + {{-109.05, 41.00}, {-109.045, 37.0}, {-102.05, 37.0}, {-102.05, 41.00}, {-109.05, 41.00}}, +} + +ctx := context.Background() +client := datasets.NewClient() + +dataset, err := client.Datasets.Get(ctx, "open_data.copernicus.sentinel2_msi") +if err != nil { + log.Fatalf("Failed to get dataset: %v", err) +} + +collection, err := client.Collections.Get(ctx, dataset.ID, "S2A_S2MSI2A") +if err != nil { + log.Fatalf("Failed to get collection: %v", err) +} + +var datapoints []*v1.Sentinel2Msi +err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, + &datapoints, + datasets.WithTemporalExtent(query.NewTimeInterval(startDate, endDate)), + datasets.WithSpatialExtent(area), +) +if err != nil { + log.Fatalf("Failed to query datapoints: %v", err) +} +``` + + +## Intersection mode + +By default, the query will return all datapoints that intersect with the specified geometry. +For certain use cases you might want to change this behavior and only return datapoints that are fully contained within the specified geometry. +Tilebox supports this behavior by specifying a mode for the spatial filter. + + + + Query results with intersects mode + Query results with intersects mode + + + Query results with contains mode + Query results with contains mode + + + +### Intersects + +The `intersects` mode is the default behavior of spatial queries. It matches all datapoints with geometries that intersect with the query geometry. + + +```python Python +area = Polygon( # area roughly covering the state of Colorado + ((-109.05, 41.00), (-109.045, 37.0), (-102.05, 37.0), (-102.05, 41.00), (-109.05, 41.00)), +) + +data = collection.query( + temporal_extent=("2025-04-02", "2025-04-03"), + # intersects is the default, so can also be omitted entirely + spatial_extent={"geometry": area, "mode": "intersects"}, +) +print(f"There are {data.sizes['time']} Sentinel-2A granules intersecting the area of Colorado on April 2nd, 2025") +``` +```go Go +startDate := time.Date(2025, 4, 2, 0, 0, 0, 0, time.UTC) +endDate := time.Date(2025, 4, 3, 0, 0, 0, 0, time.UTC) +area := orb.Polygon{ // area roughly covering the state of Colorado + {{-109.05, 41.00}, {-109.045, 37.0}, {-102.05, 37.0}, {-102.05, 41.00}, {-109.05, 41.00}}, +} + +var datapoints []*examplesv1.Sentinel2Msi +err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, + &datapoints, + datasets.WithTemporalExtent(query.NewTimeInterval(startDate, endDate)), + datasets.WithSpatialExtentFilter(&query.SpatialFilter{ + Geometry: area, + // intersects is the default, so can also be omitted entirely + Mode: datasetsv1.SpatialFilterMode_SPATIAL_FILTER_MODE_INTERSECTS, + }), +) +if err != nil { + log.Fatalf("Failed to query datapoints: %v", err) +} +``` + + +```plaintext Output +There are 27 Sentinel-2A granules intersecting the area of Colorado on April 2nd, 2025 +``` + +### Contains + + +```python Python +area = Polygon( # area roughly covering the state of Colorado + ((-109.05, 41.00), (-109.045, 37.0), (-102.05, 37.0), (-102.05, 41.00), (-109.05, 41.00)), +) + +data = collection.query( + temporal_extent=("2025-04-01", "2025-05-02"), + spatial_extent={"geometry": area, "mode": "contains"}, +) +print(f"There are {data.sizes['time']} Sentinel-2A granules fully contained within the area of Colorado on April 2nd, 2025") +``` +```go Go +startDate := time.Date(2025, 4, 2, 0, 0, 0, 0, time.UTC) +endDate := time.Date(2025, 4, 3, 0, 0, 0, 0, time.UTC) +area := orb.Polygon{ // area roughly covering the state of Colorado + {{-109.05, 41.00}, {-109.045, 37.0}, {-102.05, 37.0}, {-102.05, 41.00}, {-109.05, 41.00}}, +} + +var datapoints []*examplesv1.Sentinel2Msi +err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, + &datapoints, + datasets.WithTemporalExtent(query.NewTimeInterval(startDate, endDate)), + datasets.WithSpatialExtentFilter(&query.SpatialFilter{ + Geometry: area, + Mode: datasetsv1.SpatialFilterMode_SPATIAL_FILTER_MODE_CONTAINS, + }), +) +if err != nil { + log.Fatalf("Failed to query datapoints: %v", err) +} +``` + + +```plaintext Output +There are 16 Sentinel-2A granules fully contained within the area of Colorado on April 2nd, 2025 +``` + +## Antimeridian Crossings + +In many applications, geometries that cross the antimeridian cause issues. Since such geometries are common in satellite +data, Tilebox does take extra care to handle them out of the box correctly, by building the necessary internal spatial +index structures in a way that correctly handles antimeridian crossings and pole coverings. + +To get accurate results also at query time, it's recommend to use the `spherical` [coordinate reference system](#coordinate-reference-system) +for querying (which is the default), as it correctly handles the non-linearity introduced by the antimeridian in `cartesian` space. + + + Even if you stick to the `spherical` coordinate reference system when querying, it's still recommended to follow the + [best practices for handling geometries](/datasets/geometries). In doing so, + you can ensure that no geometry related issues will arise even when interfacing with other libraries and tools that may not properly + support non-linearities in geometries. + + +## Coordinate reference system + +Geometry intersection and containment checks can either be performed in a [3D Spherical coordinate system](https://en.wikipedia.org/wiki/Spherical_coordinate_system) +or in a standard 2D cartesian `lat/lon` coordinate system. + + + + Spherical coordinate reference system + Spherical coordinate reference system + + + Cartesian coordinate reference system + Cartesian coordinate reference system + + + +### Spherical + +The `spherical` coordinate reference system is the default and recommended one to use. It correctly handles antimeridian crossings +and is the most robust option, no matter how the datapoint geometries are cut along the antimeridian. + + + Irregardless of the coordinate reference system is used, it is always recommended to follow the best practices + for handling antimeridian crossings as described in the [Antimeridian Crossings section](#antimeridian-crossings) below. + + +When querying with the `spherical` coordinate reference system, Tilebox automatically converts all geometries to +their `x, y, z` coordinates on the unit sphere and performs the intersection and containment checks in 3D. + + +```python Python +area = Polygon( # area roughly covering the state of Colorado + ((-109.05, 41.00), (-109.045, 37.0), (-102.05, 37.0), (-102.05, 41.00), (-109.05, 41.00)), +) + +data = collection.query( + temporal_extent=("2025-04-01", "2025-05-02"), + # spherical is the default, so can also be omitted entirely + spatial_extent={"geometry": area, "coordinate_system": "spherical"}, +) +``` +```go Go +startDate := time.Date(2025, 4, 2, 0, 0, 0, 0, time.UTC) +endDate := time.Date(2025, 4, 3, 0, 0, 0, 0, time.UTC) +area := orb.Polygon{ // area roughly covering the state of Colorado + {{-109.05, 41.00}, {-109.045, 37.0}, {-102.05, 37.0}, {-102.05, 41.00}, {-109.05, 41.00}}, +} + +var datapoints []*examplesv1.Sentinel2Msi +err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, + &datapoints, + datasets.WithTemporalExtent(query.NewTimeInterval(startDate, endDate)), + datasets.WithSpatialExtentFilter(&query.SpatialFilter{ + Geometry: area, + // spherical is the default, so can also be omitted entirely + CoordinateSystem: datasetsv1.SpatialCoordinateSystem_SPATIAL_COORDINATE_SYSTEM_SPHERICAL, + }), +) +if err != nil { + log.Fatalf("Failed to query datapoints: %v", err) +} +``` + + +### Cartesian + +Tilebox can also be configured to use a standard 2D cartesian `lat/lon` coordinate system for geometry intersection and containment checks +as well. +This can be done by specifying the `cartesian` coordinate reference system when querying. + + +```python Python +area = Polygon( # area roughly covering the state of Colorado + ((-109.05, 41.00), (-109.045, 37.0), (-102.05, 37.0), (-102.05, 41.00), (-109.05, 41.00)), +) + +data = collection.query( + temporal_extent=("2025-04-01", "2025-05-02"), + spatial_extent={"geometry": area, "coordinate_system": "cartesian"}, +) +``` +```go Go +startDate := time.Date(2025, 4, 2, 0, 0, 0, 0, time.UTC) +endDate := time.Date(2025, 4, 3, 0, 0, 0, 0, time.UTC) +area := orb.Polygon{ // area roughly covering the state of Colorado + {{-109.05, 41.00}, {-109.045, 37.0}, {-102.05, 37.0}, {-102.05, 41.00}, {-109.05, 41.00}}, +} + +var datapoints []*examplesv1.Sentinel2Msi +err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, + &datapoints, + datasets.WithTemporalExtent(query.NewTimeInterval(startDate, endDate)), + datasets.WithSpatialExtentFilter(&query.SpatialFilter{ + Geometry: area, + CoordinateSystem: datasetsv1.SpatialCoordinateSystem_SPATIAL_COORDINATE_SYSTEM_CARTESIAN, + }), +) +if err != nil { + log.Fatalf("Failed to query datapoints: %v", err) +} +``` + + + + When using the `cartesian` coordinate system, antimeridian crossings may cause issues if datapoint geometries + or the query geometry do not properly respect the antimeridian cut. Check out the + [Antimeridian Crossings](#antimeridian-crossings) section below for best practices to ensure correct results irregardless + of the coordinate reference system used. + + diff --git a/datasets/query/filter-by-time.mdx b/datasets/query/filter-by-time.mdx new file mode 100644 index 0000000..d263308 --- /dev/null +++ b/datasets/query/filter-by-time.mdx @@ -0,0 +1,326 @@ +--- +title: Querying by temporal extent +sidebarTitle: Filter by time +description: Learn how to run queries for datapoints within a given temporal extent. +icon: timeline +--- + +Both [Timeseries](/datasets/types/timeseries) and [Spatio-temporal](/datasets/types/spatiotemporal) datasets support efficient time-based queries. + +To query data from a collection, use the [query](/api-reference/python/tilebox.datasets/Collection.query) method. It requires a temporal extent parameter to specify the time or time interval for querying. The behavior of the `query` method depends on the exact temporal extent parameter you provide. + +## Time interval queries + +To query data for a specific time interval, use a `tuple` in the form `(start, end)` as the `temporal_extent` parameter. Both `start` and `end` must be [TimeScalars](#time-scalar-queries), which can be `datetime` objects or strings in ISO 8601 format. + + +```python Python +from tilebox.datasets import Client + +client = Client() +sentinel2_msi = client.dataset("open_data.copernicus.sentinel2_msi") +collection = sentinel2_msi.collection("S2A_S2MSI2A") + +interval = ("2025-05-01", "2025-06-01") +data = collection.query(temporal_extent=interval, show_progress=True) +``` +```go Go +startDate := time.Date(2025, time.May, 1, 0, 0, 0, 0, time.UTC) +endDate := time.Date(2025, time.June, 1, 0, 0, 0, 0, time.UTC) +interval := query.NewTimeInterval(startDate, endDate) + +ctx := context.Background() +client := datasets.NewClient() + +dataset, err := client.Datasets.Get(ctx, "open_data.copernicus.sentinel2_msi") +if err != nil { + log.Fatalf("Failed to get dataset: %v", err) +} + +collection, err := client.Collections.Get(ctx, dataset.ID, "S2A_S2MSI2A") +if err != nil { + log.Fatalf("Failed to get collection: %v", err) +} + +var datapoints []*v1.Sentinel2Msi +err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, + &datapoints, + datasets.WithTemporalExtent(interval), +) +if err != nil { + log.Fatalf("Failed to query datapoints: %v", err) +} + +log.Printf("Queried %d datapoints", len(datapoints)) +``` + + +Output + + +```plaintext Python + Size: 39MB +Dimensions: (time: 87121) +Coordinates: + * time (time) datetime64[ns] 697kB 2025-05-01T00:00:51.02... +Data variables: (12/23) + id (time) + + + The `show_progress` parameter is optional and can be used to display a [tqdm](https://tqdm.github.io/) progress bar while loading data. + + +A time interval specified as a tuple is interpreted as a half-closed interval. This means the start time is inclusive, and the end time is exclusive. +For instance, using an end time of `2025-06-01` includes data points up to `2025-05-31 23:59:59.999`, but excludes those from `2025-06-01 00:00:00.000`. +This behavior mimics the Python `range` function and is useful for chaining time intervals. + + +```python Python +import xarray as xr + +data = [] +for month in [4, 5, 6]: + interval = (f"2025-{month}-01", f"2025-{month+1}-01") + data.append(collection.query(temporal_extent=interval, show_progress=True)) + +# Concatenate the data into a single dataset, which is equivalent +# to the result of the single request in the code example above. +data = xr.concat(data, dim="time") +``` +```go Go +var datapoints []*v1.Sentinel2Msi + +for month := 4; month <= 6; month++ { + startDate := time.Date(2025, month, 1, 0, 0, 0, 0, time.UTC) + endDate := time.Date(2025, month + 1, 1, 0, 0, 0, 0, time.UTC) + + var partialDatapoints []*v1.Sentinel2Msi + err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, + &partialDatapoints, + datasets.WithTemporalExtent(query.NewTimeInterval(startDate, endDate)), + ) + if err != nil { + log.Fatalf("Failed to query datapoints: %v", err) + } + + // Concatenate the data into a single dataset, which is equivalent + // to the result of the single request in the code example above. + datapoints = append(datapoints, partialDatapoints...) +} +``` + + +Above example demonstrates how to split a large time interval into smaller chunks while loading data in separate requests. Typically, this is not necessary as the datasets client auto-paginates large intervals. + +### Endpoint inclusivity + +For greater control over inclusivity of start and end times, you can explicitly specify a `TimeInterval`. This way you can specify both the `start` and `end` times, as well as their inclusivity. Here's an example of creating equivalent `TimeInterval` objects in two different ways. + + +```python Python +from datetime import datetime +from tilebox.datasets.data import TimeInterval + +interval1 = TimeInterval( + datetime(2021, 1, 1), datetime(2023, 1, 1), + end_inclusive=False +) +interval2 = TimeInterval( + # python datetime granularity is in milliseconds + datetime(2021, 1, 1), datetime(2022, 12, 31, 23, 59, 59, 999999), + end_inclusive=True +) + +print("Inclusivity is indicated by interval notation: ( and [") +print(interval1) +print(interval2) +print(f"They are equivalent: {interval1 == interval2}") +print(interval2.to_half_open()) + +# Query data for a time interval +data = collection.query(temporal_extent=interval1, show_progress=True) +``` +```go Go +interval1 := query.TimeInterval{ + Start: time.Date(2021, time.January, 1, 0, 0, 0, 0, time.UTC), + End: time.Date(2023, time.January, 1, 0, 0, 0, 0, time.UTC), + EndInclusive: false, +} + +interval2 := query.TimeInterval{ + Start: time.Date(2021, time.January, 1, 0, 0, 0, 0, time.UTC), + // the granularity of time.Time in Go is nanoseconds + End: time.Date(2022, time.December, 31, 23, 59, 59, 999999999, time.UTC), + EndInclusive: true, +} + +log.Println("Inclusivity is indicated by interval notation: ( and [") +log.Println(interval1.String()) +log.Println(interval2.String()) +log.Println("They are equivalent:", interval1.Equal(&interval2)) +log.Println(interval2.ToHalfOpen().String()) + +// Query data for a time interval +var datapoints []*v1.Sentinel2Msi +err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, + &datapoints, + datasets.WithTemporalExtent(interval1), +) +``` + + +```plaintext Output +Inclusivity is indicated by interval notation: ( and [ +[2021-01-01T00:00:00.000 UTC, 2023-01-01T00:00:00.000 UTC) +[2021-01-01T00:00:00.000 UTC, 2022-12-31T23:59:59.999 UTC] +They are equivalent: True +[2021-01-01T00:00:00.000 UTC, 2023-01-01T00:00:00.000 UTC) +``` + +## Time scalar queries + +You can query all datapoints linked to a specific timestamp by specifying a `TimeScalar` as the time query argument. A `TimeScalar` can be a `datetime` object or a string in ISO 8601 format. + + + Tilebox uses millisecond precision for time indexing datapoints. Thus, querying a specific time scalar, is equivalent to a time interval query of length 1 millisecond. + + +Here's how to query a data point at a specific millisecond from a [collection](/datasets/concepts/collections). + + +```python Python +data = collection.query(temporal_extent="2025-06-15T02:31:41.024") +print(data) +``` +```go Go +temporalExtent := query.NewPointInTime(time.Date(2025, time.June, 15, 2, 31, 41, 024000000, time.UTC)) + +var datapoints []*v1.Sentinel2Msi +err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, &datapoints, + datasets.WithTemporalExtent(temporalExtent), +) + +log.Printf("Queried %d datapoints", len(datapoints)) +log.Printf("First datapoint time: %s", datapoints[0].GetTime().AsTime()) +``` + + +Output + + +```plaintext Python + Size: 158kB +Dimensions: (time: 357) +Coordinates: + * time (time) datetime64[ns] 3kB 2025-06-15T02:31:41.0240... +Data variables: (12/23) + id (time) + + + A collection may contain multiple datapoints for the same millisecond, so multiple data points may be returned. If you want to fetch only a single data point, [query the collection by id](#loading-a-data-point-by-id) instead. + + + +## Timezone handling + +All `TimeScalars` specified as a string are treated as UTC if they do not include a timezone suffix. If you want to query data for a specific time or time range +in another timezone, it's recommended to a type that includes timezone information. Tilebox will automatically convert such objects to `UTC` in order to send the right query requests. +All outputs will always contain UTC timestamps, which will need to be converted again to a different timezone if required. + + +```python Python +from datetime import datetime +import pytz + +# Tokyo has a UTC+9 hours offset, so this is the same as +# 2017-01-01 02:45:25.679 UTC +tokyo_time = pytz.timezone('Asia/Tokyo').localize( + datetime(2021, 1, 1, 11, 45, 25, 679000) +) +print(tokyo_time) +data = collection.query(temporal_extent=tokyo_time) +print(data) +``` +```go Go +// Tokyo has a UTC+9 hours offset, so this is the same as +// 2017-01-01 02:45:25.679 UTC +location, _ := time.LoadLocation("Asia/Tokyo") +tokyoTime := query.NewPointInTime(time.Date(2021, 1, 1, 11, 45, 25, 679000000, location)) +log.Println(tokyoTime) + +var datapoints []*v1.Sentinel2Msi +err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, &datapoints, + datasets.WithTemporalExtent(tokyoTime), +) +if err != nil { + log.Fatalf("Failed to query datapoints: %v", err) +} + +log.Printf("Queried %d datapoints", len(datapoints)) +// time is in UTC since API always returns UTC timestamps +log.Printf("First datapoint time: %s", datapoints[0].GetTime().AsTime()) +``` + + +Output + + +```plaintext Python +2021-01-01 11:45:25.679000+09:00 + Size: 725B +Dimensions: (time: 1, latlon: 2) +Coordinates: + ingestion_time (time) datetime64[ns] 8B 2024-06-21T11:03:33.852435 + id (time) + diff --git a/datasets/query/querying-data.mdx b/datasets/query/querying-data.mdx new file mode 100644 index 0000000..66b2a2c --- /dev/null +++ b/datasets/query/querying-data.mdx @@ -0,0 +1,239 @@ +--- +title: Querying data +sidebarTitle: Querying data +description: Learn how to query and load data from Tilebox dataset collections. +icon: magnifying-glass +--- + +Tilebox offers a powerful and flexible querying API to access and filter data from your datasets. When querying, you can +[filter by time](/datasets/query/filter-by-time) and for [Spatio-temporal datasets](/datasets/types/spatiotemporal) optionally also [filter by a location in the form of a geometry](/datasets/query/filter-by-location). + +## Selecting a collection + +Querying is always done on a [collection](/datasets/concepts/collections) level, so to get started first select a collection to query. + + +```python Python +from tilebox.datasets import Client + +client = Client() +sentinel2_msi = client.dataset("open_data.copernicus.sentinel2_msi") +collection = sentinel2_msi.collection("S2A_S2MSI2A") +``` +```go Go +package main + +import ( + "context" + "log" + + "github.com/tilebox/tilebox-go/datasets/v1" +) + +func main() { + ctx := context.Background() + client := datasets.NewClient() + + dataset, err := client.Datasets.Get(ctx, "open_data.copernicus.sentinel2_msi") + if err != nil { + log.Fatalf("Failed to get dataset: %v", err) + } + + collection, err := client.Collections.Get(ctx, dataset.ID, "S2A_S2MSI2A") + if err != nil { + log.Fatalf("Failed to get collection: %v", err) + } +} +``` + + + + Querying multiple dataset collections at once is a feature already on our roadmap. If you need this functionality, please [get in touch](mailto:support@tilebox.com) so we can let you know as soon as it is available. + + +## Running a query + +To query data points from a dataset collection, use the `query` method which is available for both [python](/api-reference/python/tilebox.datasets/Collection.query) and [go](/api-reference/go/datasets/Datapoints.Query). + +Below is a simple example of querying all Sentinel-2 `S2A_S2MSI2A` data for April 2025 over the state of Colorado. + + +```python Python +from shapely import Polygon +from tilebox.datasets import Client + +area = Polygon( # area roughly covering the state of Colorado + ((-109.05, 41.00), (-109.045, 37.0), (-102.05, 37.0), (-102.05, 41.00), (-109.05, 41.00)), +) + +collection = sentinel2_msi.collection("S2A_S2MSI2A") +data = collection.query( + temporal_extent=("2025-04-01", "2025-05-01"), + spatial_extent=area, + show_progress=True, +) +``` +```go Go +startDate := time.Date(2025, 4, 1, 0, 0, 0, 0, time.UTC) +endDate := time.Date(2025, 5, 2, 0, 0, 0, 0, time.UTC) +timeInterval := query.NewTimeInterval(startDate, endDate) +area := orb.Polygon{ // area roughly covering the state of Colorado + {{-109.05, 41.00}, {-109.045, 37.0}, {-102.05, 37.0}, {-102.05, 41.00}, {-109.05, 41.00}}, +} + +collection, err := client.Collections.Get(ctx, dataset.ID, "S2A_S2MSI2A") +if err != nil { + log.Fatalf("Failed to get collection: %v", err) +} + +var datapoints []*v1.Sentinel2Msi +err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, + &datapoints, + datasets.WithTemporalExtent(timeInterval), + datasets.WithSpatialExtent(area), +) +if err != nil { + log.Fatalf("Failed to query datapoints: %v", err) +} +``` + + +To learn more about how to specify filters to narrow down the query results, check out the following sections about filtering by +time, by geometry or by datapoint ID. + + + + Learn how to run queries for datapoints within a given temporal extent. + + + Learn how to filter your query results by a certain geographical extent. + + + Find individual datapoints by their unique ID. + + + +## Automatic pagination + +Querying large datasets can result in a large number of data points. For those cases Tilebox +automatically handles pagination for you by sending paginated requests to the server. + + +When using the python SDK in an interactive notebook environment, you can additionally also display a +progress bar to keep track of the progress of the query by setting the `show_progress` parameter to `True`. + + + +## Skipping data fields + +Sometimes, only the ID or timestamp associated with a datapoint is required. In those cases, you can speed up +querying by skipping downloading of all dataset fields except of the `time`, the `id` and the `ingestion_time` by setting the `skip_data` parameter to `True`. + +For example, when checking how many datapoints exist in a given time interval, you can use `skip_data=True` to avoid loading the data fields. + + +```python Python +interval = ("2023-01-01", "2023-02-01") +data = collection.query(temporal_extent=interval, skip_data=True) +print(f"Found {data.sizes['time']} data points.") +``` +```go Go +startDate := time.Date(2023, time.January, 1, 0, 0, 0, 0, time.UTC) +endDate := time.Date(2023, time.February, 1, 0, 0, 0, 0, time.UTC) +interval := query.NewTimeInterval(startDate, endDate) + +var datapoints []*v1.Sentinel1Sar +err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, + &datapoints, + datasets.WithTemporalExtent(temporalExtent), + datasets.WithSkipData(), +) +if err != nil { + log.Fatalf("Failed to query datapoints: %v", err) +} + +log.Printf("Queried %d datapoints", len(datapoints)) +log.Printf("First datapoint time: %s", datapoints[0].GetTime().AsTime()) +``` + + +Output + + +```plaintext Python + Size: 160B +Dimensions: (time: 1) +Coordinates: + ingestion_time (time) datetime64[ns] 8B 2024-08-01T08:53:08.450499 + id (time) + +## Empty response + +Query will not raise an error if no data points are found for the specified query, instead an empty result is returned. + +In python, this is an empty `xarray.Dataset` object. In Go, an empty slice of datapoints. To check for an empty response, you can coerce the result to a boolean +or check the length of the slice. + + +```python Python +timestamp_with_no_data_points = "1997-02-06 10:21:00:000" +data = collection.query(temporal_extent=timestamp_with_no_data_points) +if not data: + print("No data points found") +``` +```go Go +timeWithNoDatapoints := query.NewPointInTime(time.Date(1997, time.February, 6, 10, 21, 0, 0, time.UTC)) + +var datapoints []*v1.Sentinel1Sar +err = client.Datapoints.QueryInto(ctx, + []uuid.UUID{collection.ID}, + &datapoints, + datasets.WithTemporalExtent(timeWithNoDatapoints), +) +if err != nil { + log.Fatalf("Failed to query datapoints: %v", err) +} + +if len(datapoints) == 0 { + log.Println("No data points found") +} +``` + + +Output + + +```plaintext Python +No data points found +``` +```plaintext Go +No data points found +``` + + diff --git a/datasets/types/spatiotemporal.mdx b/datasets/types/spatiotemporal.mdx index 10acb1f..25340bd 100644 --- a/datasets/types/spatiotemporal.mdx +++ b/datasets/types/spatiotemporal.mdx @@ -4,10 +4,6 @@ description: Spatio-temporal datasets link each data point to a specific point i icon: earth-europe --- - - Spatio-temporal datasets are currently in development and not available yet. Stay tuned for updates - - Each spatio-temporal dataset comes with a set of required and auto-generated fields for each data point. ## Required fields @@ -23,7 +19,7 @@ While the specific data fields between different time series datasets can vary, - A location on the earth's surface associated with each data point. Supported geometry types are `Polygon`, `MultiPolygon`, `Point` and `MultiPoint`. + A location on the earth's surface associated with each data point. Supported geometry types are `Point`, `Polygon` and `MultiPolygon`. ## Auto-generated fields @@ -48,4 +44,9 @@ already outlined will be automatically added to the dataset schema. ## Spatio-temporal queries Spatio-temporal datasets support efficient time-based and spatially filtered queries. To query a specific location in a given time interval, -specify a time range and a geometry when [querying data points](/datasets/query) from a collection. +specify a time range and a geometry when [querying data points](/datasets/query/filter-by-location) from a collection. + +## Geometries + +Handling Geometries can traditionally be a bit tricky, especially when working with geometries that cross the antimeridian or cover a pole. +Tilebox is designed to take away most of the friction involved in this, but it's still recommended to follow the [best practices for handling geometries](/datasets/geometries). diff --git a/datasets/types/timeseries.mdx b/datasets/types/timeseries.mdx index 77c2034..5d05822 100644 --- a/datasets/types/timeseries.mdx +++ b/datasets/types/timeseries.mdx @@ -39,4 +39,4 @@ already outlined will be automatically added to the dataset schema. ## Time-based queries -Timeseries datasets support time-based queries. To query a specific time interval, specify a time range when [querying data](/datasets/query) from a collection. +Timeseries datasets support time-based queries. To query a specific time interval, specify a time range when [querying data](/datasets/query/filter-by-time) from a collection. diff --git a/docs.json b/docs.json index 0ee7fa5..6ead708 100644 --- a/docs.json +++ b/docs.json @@ -43,9 +43,19 @@ "datasets/types/spatiotemporal" ] }, - "datasets/query", + { + "group": "Query", + "icon": "server", + "pages": [ + "datasets/query/querying-data", + "datasets/query/filter-by-time", + "datasets/query/filter-by-location", + "datasets/query/filter-by-id" + ] + }, "datasets/ingest", "datasets/delete", + "datasets/geometries", "datasets/open-data" ] }, @@ -126,8 +136,7 @@ "sdks/python/install", "sdks/python/sample-notebooks", "sdks/python/xarray", - "sdks/python/async", - "sdks/python/geometries" + "sdks/python/async" ] }, { diff --git a/guides/datasets/create.mdx b/guides/datasets/create.mdx index 25a4e90..54c4496 100644 --- a/guides/datasets/create.mdx +++ b/guides/datasets/create.mdx @@ -1,6 +1,6 @@ --- title: Creating a dataset -description: Learn how to create a dataset +description: Learn how to create a dataset with your own, custom schema icon: database --- @@ -47,8 +47,9 @@ This page guides you through the process of creating a dataset in Tilebox using Specify the fields for your dataset. Each field has these properties: - - `Name` is the name of the field (it should be `snake_case`). - - `Type` and `Array` let you specify the field data type and whether it's an array. See below for an explanation of the available data. + - `Name` is the name of the field (by convention it is recommended to be in `snake_case`). + - `Type` is the data type of the field. + - `Array` can be set to indicate that the field contains multiple values of the specified type. - `Description` is an optional brief description of the field. You can use it to provide more context and details about the data. - `Example value` is an optional example for this field. It can be useful for documentation purposes. @@ -96,5 +97,5 @@ Once you are done editing the documentation, click the `Save` button to save you You can always add new fields to a dataset. -If you want to remove or edit existing fields, you'll first need to empty all collections in the dataset. -Then, you can edit the dataset schema in the console. +If you want to remove or edit existing fields, you'll first need to empty all collections in the dataset, to ensure that +no existing data points relying on those fields exist. Then, you can freely edit the dataset schema in the console. diff --git a/guides/datasets/ingest-format.mdx b/guides/datasets/ingest-format.mdx index 47d824b..1119ad5 100644 --- a/guides/datasets/ingest-format.mdx +++ b/guides/datasets/ingest-format.mdx @@ -4,10 +4,18 @@ description: Learn how to ingest data from common file formats into Tilebox icon: file-binary --- -Through the usage of `xarray` and `pandas` you can also easily ingest existing datasets available in file -formats, such as CSV, [Parquet](https://parquet.apache.org/), [Feather](https://arrow.apache.org/docs/python/feather.html) and more. +For ingesting data from common file formats, it's recommend to use the [Tilebox Python SDK](/sdks/python/install), since +it provides out-of-the-box support for reading many common formats through third party libraries for loading data as either +`pandas.DataFrame` or `xarray.Dataset`, which can then be directly ingested into Tilebox. -## CSV + +## Reading and previewing the data + +To ingest data from a file, you first need to read it into a `pandas.DataFrame` or an `xarray.Dataset`. +How that can be achieved depends on the file format. The following sections show examples for a couple of common +file formats. + +### CSV Comma-separated values (CSV) is a common file format for tabular data. It's widely used in data science. Tilebox supports CSV ingestion using the `pandas.read_csv` function. @@ -21,22 +29,15 @@ time,value,sensor,precise_time,sensor_history,some_unwanted_column 2025-03-28T11:45:19Z,273.15,B,2025-03-28T11:45:19.128742312Z,"[300.16, 280.12, 273.15]","Unsupported" ``` -This data already conforms to the schema of the `MyCustomDataset` dataset, except for `some_unwanted_column` which -you want to drop before you ingest it. Here is how this could look like: - ```python Python import pandas as pd data = pd.read_csv("ingestion_data.csv") -data = data.drop(columns=["some_unwanted_column"]) - -collection = dataset.get_or_create_collection("CSVMeasurements") -collection.ingest(data) ``` -## Parquet +### Parquet [Apache Parquet](https://parquet.apache.org/) is an open source, column-oriented data file format designed for efficient data storage and retrieval. Tilebox supports Parquet ingestion using the `pandas.read_parquet` function. @@ -48,19 +49,32 @@ The parquet file used in this example [is available here](https://storage.google import pandas as pd data = pd.read_parquet("ingestion_data.parquet") +``` + -# our data already conforms to the schema of the MyCustomDataset -# dataset, so lets ingest it -collection = dataset.get_or_create_collection("ParquetMeasurements") -collection.ingest(data) +### GeoParquet + +[GeoParquet](https://geoparquet.org/) is an extension of the Parquet file format, adding geospatial +features support to Parquet. Tilebox supports GeoParquet ingestion using the `geopandas.read_parquet` function. + +The GeoParquet file used in this example [is available here](https://storage.googleapis.com/tbx-web-assets-2bad228/docs/data-samples/modis_MCD12Q1.geoparquet). + + +```python Python +import geopandas as gpd + +data = gpd.read_parquet("modis_MCD12Q1.geoparquet") ``` -## Feather + + For a step-by-step guide of ingesting a GeoParquet file, check out our [Ingesting data](/guides/datasets/ingest) guide. + + +### Feather [Feather](https://arrow.apache.org/docs/python/feather.html) is a file format originating from the Apache Arrow project, -designed for storing tabular data in a fast and memory-efficient way. It's supported by many programming languages, -including Python. Tilebox supports Feather ingestion using the `pandas.read_feather` function. +designed for storing tabular data in a fast and memory-efficient way. Tilebox supports Feather ingestion using the `pandas.read_feather` function. The feather file used in this example [is available here](https://storage.googleapis.com/tbx-web-assets-2bad228/docs/data-samples/ingestion_data.feather). @@ -69,14 +83,56 @@ The feather file used in this example [is available here](https://storage.google import pandas as pd data = pd.read_feather("ingestion_data.feather") +``` + -# our data already conforms to the schema of the MyCustomDataset -# dataset, so lets ingest it -collection = dataset.get_or_create_collection("FeatherMeasurements") -collection.ingest(data) +## Mapping columns to dataset fields + +Once data is read into a `pandas.DataFrame` or an `xarray.Dataset`, it can be ingested into Tilebox directly. +The column names of the `pandas.DataFrame` or the variables and coordinates of the `xarray.Dataset` are mapped to the fields of the +Tilebox dataset to ingest into. + +Depending on how closely the column names or variable/coordinate names match the field names in the Tilebox dataset, +you might need to rename some columns/variables/coordinates before ingestion. + +### Renaming fields + + +```python Pandas +data = data.rename({"precise_time": "measurement_time"}) +``` +```python Xarray +data = data.rename({"precise_time": "measurement_time"}) ``` -## GeoParquet +## Dropping fields + +In case you want to skip certain columns/variables/coordinates entirely, you can drop them before ingestion. + + +```python Pandas +data = data.drop(columns=["some_unwanted_column"]) +``` +```python Xarray +data = data.drop_vars(["some_unwanted_variable"]) +``` + + +## Ingesting the data + +Once the data is in the correct format, you can ingest it into Tilebox. + + +```python Python +from tilebox.datasets import Client + +client = Client() +# my_custom_dataset has a schema that matches the data we want to ingest +dataset = client.dataset("my_org.my_custom_dataset") + +collection = dataset.get_or_create_collection("Measurements") +collection.ingest(data) +``` + -Please check out the [Ingesting data](/guides/datasets/ingest) guide for an example of ingesting a GeoParquet file. diff --git a/guides/datasets/ingest.mdx b/guides/datasets/ingest.mdx index 01fb045..a7de461 100644 --- a/guides/datasets/ingest.mdx +++ b/guides/datasets/ingest.mdx @@ -41,16 +41,16 @@ required packages using your preferred package manager. For new projects, Tilebo ```bash uv -uv add tilebox-datasets geopandas folium matplotlib mapclassify +uv add tilebox-datasets geopandas lonboard ``` ```bash pip -pip install tilebox-datasets geopandas folium matplotlib mapclassify +pip install tilebox-datasets geopandas lonboard ``` ```bash poetry -poetry add tilebox-datasets="*" geopandas="*" folium="*" matplotlib="*" mapclassify="*" +poetry add tilebox-datasets="*" geopandas="*" lonboard="*" ``` ```bash pipenv -pipenv install tilebox-datasets geopandas folium matplotlib mapclassify +pipenv install tilebox-datasets geopandas lonboard ``` @@ -84,7 +84,9 @@ Geopandas comes with a built in explorer to visually explore the dataset. ```python Python -modis_data.head(1000).explore(width=800, height=600) +from lonboard import viz + +viz(modis_data, map_kwargs={"show_tooltip": True}) ``` @@ -95,9 +97,9 @@ modis_data.head(1000).explore(width=800, height=600) ## Create a Tilebox dataset -Now you'll create a [Timeseries](/datasets/types/timeseries) dataset with the same schema as the given MODIS dataset. +Now you'll create a [Spatio-temporal](/datasets/types/spatio-temporal) dataset with the same schema as the given MODIS dataset. To do so, you'll use the [Tilebox Console](/console), navigate to `My Datasets` and click `Create Dataset`. Then select -`Timeseries Dataset` as the dataset type. +`Spatio-temporal Dataset` as the dataset type. For more information on creating a dataset, check out the [Creating a dataset](/guides/datasets/create) guide for a @@ -109,7 +111,6 @@ Now, to match the given MODIS dataset, you'll specify the following fields: | Field | Type | Note | | --- | --- | --- | | `granule_name` | string | MODIS granule name | -| `geometry` | Geometry | Tile boundary coordinates of the granule | | `end_time` | Timestamp | Measurement end time | | `horizontal_tile_number` | int64 | Horizontal modis tile number (0-35) | | `vertical_tile_number` | int64 | Vertical modis tile number (0-17) | @@ -135,8 +136,8 @@ which was assigned automatically based on the specified `code_name`. To find out in the console. - Explore the MODIS dataset - Explore the MODIS dataset + Dataset slug display in the console + Dataset slug display in the console You can now instantiate the dataset client and access the dataset. @@ -187,36 +188,45 @@ You can now query the newly ingested data. You can query a subset of the data fo ```python Python -data = collection.query(temporal_extent=("2015-01-01", "2020-01-01")) +from shapely import Polygon + +area = Polygon( # area roughly covering the US + ((-124.45, 49.19), (-120.88, 29.31), (-66.87, 24.77), (-65.34, 47.84), (-124.45, 49.19)), +) + +data = collection.query( + temporal_extent=("2015-01-01", "2020-01-01"), + spatial_extent=area +) data ``` ```plaintext Output - Size: 403kB -Dimensions: (time: 1575) + Size: 28kB +Dimensions: (time: 110) Coordinates: - * time (time) datetime64[ns] 13kB 2015-01-01 ... 2019-01-01 + * time (time) datetime64[ns] 880B 2015-01-01 ... 2019-01-01 Data variables: (12/14) - id (time) - For more information on accessing and querying data, check out [querying data](/datasets/query). + For more information on accessing and querying data, check out [querying data](/datasets/query/querying-data). ## View the data in the console diff --git a/quickstart.mdx b/quickstart.mdx index b8cc36a..1010771 100644 --- a/quickstart.mdx +++ b/quickstart.mdx @@ -192,7 +192,7 @@ If you prefer to work locally, follow these steps to get started. // load data from a collection in a given time range and spatial extent colorado := orb.Polygon{ - {{-109.05, 37.09}, {-102.06, 37.09}, {-102.06, 41.59}, {-109.05, 41.59}, {-109.05, 37.09}}, + {{-109.05, 41.00}, {-109.045, 37.0}, {-102.05, 37.0}, {-102.05, 41.00}, {-109.05, 41.00}}, } startDate := time.Date(2025, time.March, 1, 0, 0, 0, 0, time.UTC) endDate := time.Date(2025, time.April, 1, 0, 0, 0, 0, time.UTC) diff --git a/sdks/python/geometries.mdx b/sdks/python/geometries.mdx deleted file mode 100644 index 39a7533..0000000 --- a/sdks/python/geometries.mdx +++ /dev/null @@ -1,321 +0,0 @@ ---- -title: Geometries -description: How geometries are handled in the Tilebox Python client. -icon: earth-americas ---- - -Many datasets consist of granules that represent specific geographical areas on the Earth's surface. Often, a polygon defining the outline of these areas—a footprint—accompanies other granule metadata in time series datasets. Tilebox provides built-in support for working with geometries. - -Here's an example that loads some granules from the `ERS SAR` Opendata dataset, which contains geometries. - -```python Loading ERS data -from tilebox.datasets import Client - -client = Client() -datasets = client.datasets() - -ers_collection = datasets.open_data.asf.ers_sar.collection("ERS-2") -ers_data = ers_collection.query(temporal_extent=("2008-02-10T21:00", "2008-02-10T22:00")) -``` - -## Shapely - -In the `ers_data` dataset, each granule includes a `geometry` field that represents the footprint of each granule as a polygon. Tilebox automatically converts geometry fields to `Polygon` or `MultiPolygon` objects from the [Shapely](https://shapely.readthedocs.io/en/stable/manual.html) library. By integrating with Shapely, you can use the rich set of libraries and tools it provides. That includes support for computing polygon characteristics such as total area, intersection checks, and conversion to other formats. - -```python Printing geometries -geometries = ers_data.geometry.values -print(geometries) -``` - -Each geometry is a [shapely.Geometry](https://shapely.readthedocs.io/en/stable/geometry.html#geometry). - -```plaintext Output -[ - - - - - - - - - - - - ] -``` - - - Geometries are not always [Polygon](https://shapely.readthedocs.io/en/stable/reference/shapely.Polygon.html#shapely.Polygon) objects. More complex footprint geometries are represented as [MultiPolygon](https://shapely.readthedocs.io/en/stable/reference/shapely.MultiPolygon.html#shapely.MultiPolygon) objects. - - -### Accessing Coordinates - -You can select a polygon from the geometries and access the underlying coordinates and an automatically computed centroid point. - -```python Accessing coordinates and computing a centroid point -polygon = geometries[0] -lon, lat = polygon.exterior.coords.xy -center, = list(polygon.centroid.coords) - -print(lon) -print(lat) -print(center) -``` - -```plaintext Output -array('d', [-150.753244, -152.031574, -149.183655, -147.769339, -150.753244]) -array('d', [74.250081, 73.336051, 73.001748, 73.899483, 74.250081]) -(-149.92927907414239, 73.62538063474753) -``` - - - Interactive environments such as [Jupyter Notebooks](/sdks/python/sample-notebooks) can visualize Polygon shapes graphically. Just type `polygon` in an empty cell and execute it for a visual representation of the polygon shape. - - -### Visualization on a Map - -To visualize polygons on a map, you can use [Folium](https://pypi.org/project/folium/). Below is a helper function that produces an OpenStreetMap with the Polygon overlaid. - -```python visualize helper function -# pip install folium -from folium import Figure, Map, Polygon as FoliumPolygon, GeoJson, TileLayer -from folium.plugins import MiniMap -from shapely import Polygon, to_geojson -from collections.abc import Iterable - -def visualize(poly: Polygon | Iterable[Polygon], zoom=4): - """Visualize a polygon or a list of polygons on a map""" - if not isinstance(poly, Iterable): - poly = [poly] - - fig = Figure(width=600, height=600) - map = Map(location=geometries[len(geometries)//2].centroid.coords[0][::-1], zoom_start=zoom, control_scale=True) - map.add_child(MiniMap()) - fig.add_child(map) - - for p in poly: - map.add_child(GeoJson(to_geojson(p))) - return fig -``` - -Here's how to use it. - -```python Visualizing a polygon -visualize(polygon) -``` - -Single visualization - -The `visualize` helper function supports a list of polygons, which can display the data layout of the ERS granules. - -```python Visualizing multiple polygons -visualize(geometries) -``` - -Granules visualization - -## Format conversion - -Shapely supports converting Polygons to some common formats, such as [GeoJSON](https://geojson.org/) or [Well-Known Text (WKT)](https://docs.ogc.org/is/18-010r7/18-010r7.html). - -```python Converting to GeoJSON -from shapely import to_geojson - -print(to_geojson(polygon)) -``` - -```plaintext Output -{"type":"Polygon","coordinates":[[[-150.753244,74.250081],[-152.031574,73.336051],[-149.183655,73.001748],[-147.769339,73.899483],[-150.753244,74.250081]]]} -``` - -```python Converting to WKT -from shapely import to_wkt - -print(to_wkt(polygon)) -``` - -```plaintext Output -POLYGON ((-150.753244 74.250081, -152.031574 73.336051, -149.183655 73.001748, -147.769339 73.899483, -150.753244 74.250081)) -``` - -## Checking intersections - -One common task when working with geometries is checking if a given geometry falls into a specific area of interest. Shapely provides an `intersects` method for this purpose. - -```python Checking intersections -from shapely import box - -# Box representing the rectangular area lon=(-160, -150) and lat=(69, 70) -area_of_interest = box(-160, 69, -150, 70) - -for i, polygon in enumerate(geometries): - if area_of_interest.intersects(polygon): - print(f"{ers_data.granule_name[i].item()} intersects the area of interest!") - else: - print(f"{ers_data.granule_name[i].item()} doesn't intersect the area of interest!") -``` - -```plaintext Output -E2_66974_STD_F264 doesn't intersect the area of interest! -E2_66974_STD_F265 doesn't intersect the area of interest! -E2_66974_STD_F267 doesn't intersect the area of interest! -E2_66974_STD_F269 doesn't intersect the area of interest! -E2_66974_STD_F271 doesn't intersect the area of interest! -E2_66974_STD_F273 intersects the area of interest! -E2_66974_STD_F275 intersects the area of interest! -E2_66974_STD_F277 intersects the area of interest! -E2_66974_STD_F279 doesn't intersect the area of interest! -E2_66974_STD_F281 doesn't intersect the area of interest! -E2_66974_STD_F283 doesn't intersect the area of interest! -E2_66974_STD_F285 doesn't intersect the area of interest! -E2_66974_STD_F289 doesn't intersect the area of interest! -``` - -## Combining polygons - -As shown in the visualization of the granule footprints, the granules collectively form an orbit from pole to pole. Measurements are often combined during processing. You can do the same with geometries by combining them into a single polygon, which represents the hull around all individual footprints using [shapely.unary_union](https://shapely.readthedocs.io/en/stable/reference/shapely.unary_union.html). - -```python Combining multiple polygons -from shapely.ops import unary_union - -hull = unary_union(geometries) -visualize(hull) -``` - -The computed hull consists of two polygons due to a gap (probably a missing granule) in the geometries. Such geometries are represented as [Multi Polygons](#multi-polygons). - -Unary union visualization - -## Multi Polygons - -A collection of one or more non-overlapping polygons combined into a single geometry is called a [MultiPolygon](https://shapely.readthedocs.io/en/latest/reference/shapely.MultiPolygon.html). Footprint geometries can be of type `MultiPolygon` due to gaps or pole discontinuities. The computed hull in the previous example is a `MultiPolygon`. - -```python Accessing individual polygons of a MultiPolygon -print(f"The computed hull of type {type(hull).__name__} consists of {len(hull.geoms)} sub polygons") -for i, poly in enumerate(hull.geoms): - print(f"Sub polygon {i} has an area of {poly.area}") -``` - -```plaintext Output -The computed hull of type MultiPolygon consists of 2 sub polygons -Sub polygon 0 has an area of 2.025230449898011 -Sub polygon 1 has an area of 24.389998081651527 -``` - -## Antimeridian Crossings - -A common issue with `longitude / latitude` geometries is crossings of the 180-degree meridian, or the antimeridian. For example, the coordinates of a `LineString` from Japan to the United States might look like this: - -`140, 141, 142, ..., 179, 180, -179, -178, ..., -125, -124` - -Libraries like Shapely are not designed to handle spherical coordinate systems, so caution is necessary with such geometries. - -Here's an `ERS` granule demonstrating this issue. - -```python Antimeridian Crossing -# A granule that crosses the antimeridian -granule = ers_collection.find("0119bb86-0260-5819-6aab-f99796417155") -polygon = granule.geometry.item() -print(polygon.exterior.coords.xy) -visualize(polygon) -``` - -```plaintext Output -array('d', [177.993407, 176.605009, 179.563047, -178.904076, 177.993407]) -array('d', [74.983185, 74.074615, 73.727752, 74.61847, 74.983185]) -``` - -Antimeridian buggy visualization - -This 2D visualization appears incorrect. Both the visualization and any calculations performed may yield inaccurate results. For instance, testing whether the granule intersects the 0-degree meridian provides a false positive. - -```python Problems with calculating intersections -from shapely import LineString - -null_meridian = LineString([(0, -90), (0, 90)]) -print(polygon.intersects(null_meridian)) # True - but this is incorrect! -``` - -The GeoJSON specification offers a solution for this problem. In the section [Antimeridian Cutting](https://datatracker.ietf.org/doc/html/rfc7946#section-3.1.9), it suggests always cutting lines and polygons into two parts—one for the eastern hemisphere and one for the western hemisphere. -In python, this can be achieved using the [antimeridian](https://pypi.org/project/antimeridian/) package. - -```python Cutting the polygon along the antimeridian -# pip install antimeridian - -import antimeridian -fixed_polygon = antimeridian.fix_polygon(polygon) -visualize(fixed_polygon) - -print(fixed_polygon.intersects(null_meridian)) # False - this is correct now -``` - -Antimeridian fixed visualization - -Since Shapely is unaware of the spherical nature of this data, the **centroid** of the fixed polygon **is still incorrect**. The antimeridian package also includes a function to correct this. - -```python Calculating the centroid of a cut polygon crossing the antimeridian -print("Wrongly computed centroid coordinates (Shapely)") -print(list(fixed_polygon.centroid.coords)) -print("Correct centroid coordinates (Antimeridian taken into account)") -print(list(antimeridian.centroid(fixed_polygon).coords)) -``` - -```plaintext Output -Wrongly computed centroid coordinates (shapely) -[(139.8766350146937, 74.3747116658462)] -Correct centroid coordinates (antimeridian taken into account) -[(178.7782777050171, 74.3747116658462)] -``` - -## Spherical Geometry - -Another approach to handle the antimeridian issue is performing all coordinate-related calculations, such as polygon intersections, in a [spherical coordinate system](https://en.wikipedia.org/wiki/Spherical_coordinate_system). - -One useful library for this is [spherical_geometry](https://spherical-geometry.readthedocs.io/en/latest/). Here's an example. - -```python Spherical Geometry -# pip install spherical-geometry - -from spherical_geometry.polygon import SphericalPolygon -from spherical_geometry.vector import lonlat_to_vector - -lon, lat = polygon.exterior.coords.xy -spherical_poly = SphericalPolygon.from_lonlat(lon, lat) -# Let's check the x, y, z coordinates of the spherical polygon: -print(list(spherical_poly.points)) -``` - -```plaintext Output -[array([[-0.25894363, 0.00907234, 0.96584983], - [-0.2651968 , -0.00507317, 0.96418096], - [-0.28019363, 0.00213687, 0.95994112], - [-0.27390375, 0.01624885, 0.96161984], - [-0.25894363, 0.00907234, 0.96584983]])] -``` - -Now, you can compute intersections or check if a particular point is within the polygon. You can compare the incorrect calculation using `shapely` with the correct version when using `spherical_geometry`. - -```python Correct calculations using spherical geometry -# A point on the null-meridian, way off from our polygon -null_meridian_point = 0, 74.4 -# A point actually inside our polygon -point_inside = 178.8, 74.4 - -print("Shapely results:") -print("- Null meridian point inside:", polygon.contains(Point(*null_meridian_point))) -print("- Actual inside point inside:", polygon.contains(Point(*point_inside))) - -print("Spherical geometry results:") -print("- Null meridian point inside:", spherical_poly.contains_lonlat(*null_meridian_point)) -print("- Actual inside point inside:", spherical_poly.contains_lonlat(*point_inside)) -``` - -```plaintext Output -Shapely results: -- Null meridian point inside: True -- Actual inside point inside: False -Spherical geometry results: -- Null meridian point inside: False -- Actual inside point inside: True -```