Xarray reader

[vincentsarago]

ping @geospatial-jeff @kylebarron @TomAugspurger @rabernat

[rabernat]

🙌

[rabernat]

I recommend using this dataset: https://pangeo-forge.org/dashboard/feedstock/3

import xarray as xr

store = 'https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/noaa-coastwatch-geopolar-sst-feedstock/noaa-coastwatch-geopolar-sst.zarr'
ds = xr.open_dataset(store, engine='zarr', chunks={})
ds

[geospatial-jeff]

This seems to work:

import mercantile
import rioxarray
import xarray as xr
import numpy as np


def clip_to_tile(tile: mercantile.Tile, ds: xr.DataArray):
    bounds = mercantile.xy_bounds(tile)
    wgs84_bounds = mercantile.bounds(tile)
    
    # Create source array by clipping the xarray dataset to extent of the tile.
    source_arr = ds.rio.clip_box(*wgs84_bounds)
    
    # Create the target array.
    x = np.linspace(bounds[0], bounds[2], 256)
    y = np.linspace(bounds[1], bounds[3], 256)
    target_arr = xr.DataArray(
        np.zeros((ds['time'].shape[0], 256, 256), ds.dtype),
        [
            ("time", ds['time'].data),
            ("y", y),
            ("x", x)
        ]
    )
    target_arr.rio.write_crs(3857, inplace=True)
    
    # Map source to target
    matched = source_arr.rio.reproject_match(target_arr)
    return matched

# Create a web mercator tile.
x = 0
y = 0
z = 10
tile = mercantile.tile(x, y, z)

# Open xarray dataset.
store = 'https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/noaa-coastwatch-geopolar-sst-feedstock/noaa-coastwatch-geopolar-sst.zarr'
ds = xr.open_dataset(store, engine='zarr', chunks={})
data_array = ds['analysed_sst']

# Add CRS so rioxarray doens't complain.
data_array.rio.write_crs(4326, inplace=True)

# Only fetch one time slice.
data_array = data_array[:1]

arr = clip_to_tile(tile, data_array)

[rabernat]

Here's how you would make a cached zarr store

import zarr
url = 'https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr'
store = zarr.storage.FSStore(url)
cache = zarr.LRUStoreCache(store, max_size=2**28)
ds = xr.open_dataset(cache, engine='zarr')

[vincentsarago]

🥳

import rioxarray 
import xarray
from matplotlib.pyplot import imshow
from rio_tiler.io import XarrayReader

with xarray.open_dataset(
    "https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/noaa-coastwatch-geopolar-sst-feedstock/noaa-coastwatch-geopolar-sst.zarr",
    engine="zarr",
    decode_coords="all"
) as src:
    with XarrayReader(None, dataset=src["analysed_sst"][:1], crs="epsg:4326") as dst:
        img = dst.tile(1, 1, 2)
        imshow(img.data_as_image())

[vincentsarago]

with the latest commit

with xarray.open_dataset(
    "https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/noaa-coastwatch-geopolar-sst-feedstock/noaa-coastwatch-geopolar-sst.zarr",
    engine="zarr",
    decode_coords="all"
) as src:
    ds = src["analysed_sst"][:1]
    ds.rio.write_crs("epsg:4326", inplace=True)
    with XarrayReader(ds) as dst:
        img = dst.tile(1, 1, 2)

[TomAugspurger]

An example using Daymet on Azure:

import rioxarray 
import xarray
from matplotlib.pyplot import imshow
from rio_tiler.io import XarrayReader
import requests

credential = requests.get("https://planetarycomputer.microsoft.com/api/sas/v1/token/daymet-daily-pr").json()["token"]

url = "az://daymet-zarr/annual/pr.zarr"
storage_options = {
    "account_name": "daymeteuwest",
    "credential": credential
}
src = xarray.open_dataset(url, engine="zarr", decode_coords="all", storage_options=storage_options)

with XarrayReader(src["prcp"][:1]) as dst:
    img = dst.tile(1, 1, 2)
    imshow(img.data_as_image())

With decode_coords="all", rioxarray is able to correctly construct a CRS (I think using CF-conventions and the grid_mapping_name attribute).

[vincentsarago]

with latest commit I've removed XarrayReader from rio_tiler.io.__init__ so you need to do from rio_tiler.io.xarray import XarrayReader

[vincentsarago]

here is a simple tile server based on XarrayReader

import urllib
from enum import Enum
from typing import List, Optional, Tuple

import xarray
from fastapi import FastAPI, Query
from fastapi.responses import Response
from pydantic import BaseModel, Field, root_validator
from rio_tiler.io.xarray import XarrayReader
from starlette.requests import Request


class SchemeEnum(str, Enum):
    """TileJSON scheme choice."""

    xyz = "xyz"
    tms = "tms"


class TileJSON(BaseModel):
    """
    TileJSON model.

    Based on https://github.com/mapbox/tilejson-spec/tree/master/2.2.0

    """

    tilejson: str = "2.2.0"
    name: Optional[str]
    description: Optional[str]
    version: str = "1.0.0"
    attribution: Optional[str]
    template: Optional[str]
    legend: Optional[str]
    scheme: SchemeEnum = SchemeEnum.xyz
    tiles: List[str]
    grids: Optional[List[str]]
    data: Optional[List[str]]
    minzoom: int = Field(0, ge=0, le=30)
    maxzoom: int = Field(30, ge=0, le=30)
    bounds: List[float] = [-180, -90, 180, 90]
    center: Optional[Tuple[float, float, int]]

    @root_validator
    def compute_center(cls, values):
        """Compute center if it does not exist."""
        bounds = values["bounds"]
        if not values.get("center"):
            values["center"] = (
                (bounds[0] + bounds[2]) / 2,
                (bounds[1] + bounds[3]) / 2,
                values["minzoom"],
            )
        return values

    class Config:
        """TileJSON model configuration."""

        use_enum_values = True


app = FastAPI()


@app.get("/tiles/{z}/{x}/{y}", response_class=Response)
def tile(
    z: int,
    x: int,
    y: int,
    url: str = Query(description="Zarr URL"),
    variable: str = Query(description="Zarr Variable"),
):
    with xarray.open_dataset(url, engine="zarr", decode_coords="all") as src:
        ds = src[variable][:1]

        # Make sure we are a CRS
        crs = ds.rio.crs or "epsg:4326"
        ds.rio.write_crs(crs, inplace=True)

        with XarrayReader(ds) as dst:
            img = dst.tile(x, y, z)

        content = img.render(img_format="PNG")

    return Response(content, media_type="image/png")


@app.get(
    "/tilejson.json",
    response_model=TileJSON,
    responses={200: {"description": "Return a tilejson"}},
    response_model_exclude_none=True,
    tags=["API"],
)
def tilejson(
    request: Request,
    url: str = Query(description="Zarr URL"),
    variable: str = Query(description="Zarr Variable"),
):
    """Handle /tilejson.json requests."""
    kwargs: Dict[str, Any] = {"z": "{z}", "x": "{x}", "y": "{y}"}
    tile_url = request.url_for("tile", **kwargs)

    qs = [
        (key, value)
        for (key, value) in request.query_params._list
        if key not in ["tile_format"]
    ]
    if qs:
        tile_url += f"?{urllib.parse.urlencode(qs)}"

    with xarray.open_dataset(url, engine="zarr", decode_coords="all") as src:
        ds = src[variable][:1]

        # Make sure we are a CRS
        crs = ds.rio.crs or "epsg:4326"
        ds.rio.write_crs(crs, inplace=True)

        with XarrayReader(ds) as dst:
            return dict(
                bounds=dst.geographic_bounds,
                minzoom=dst.minzoom,
                maxzoom=dst.maxzoom,
                name="xarray",
                tilejson="2.1.0",
                tiles=[tile_url],
            )

curl http://127.0.0.1:8081/tilejson.json?url=https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/noaa-coastwatch-geopolar-sst-feedstock/noaa-coastwatch-geopolar-sst.zarr&variable=analysed_sst | jq

{
  "tilejson": "2.1.0",
  "name": "xarray",
  "version": "1.0.0",
  "scheme": "xyz",
  "tiles": [
    "[http://127.0.0.1:8081/tiles/{z}/{x}/{y}?url=https%3A%2F%2Fncsa.osn.xsede.org%2FPangeo%2Fpangeo-forge%2Fnoaa-coastwatch-geopolar-sst-feedstock%2Fnoaa-coastwatch-geopolar-sst.zarr&variable=analysed_sst](http://127.0.0.1:8081/tiles/%7Bz%7D/%7Bx%7D/%7By%7D?url=https%3A%2F%2Fncsa.osn.xsede.org%2FPangeo%2Fpangeo-forge%2Fnoaa-coastwatch-geopolar-sst-feedstock%2Fnoaa-coastwatch-geopolar-sst.zarr&variable=analysed_sst)"
  ],
  "minzoom": 0,
  "maxzoom": 2,
  "bounds": [
    -180.00000610436345,
    -89.99999847369712,
    180.00000610436345,
    89.99999847369712
  ],
  "center": [
    0,
    0,
    0
  ]
}

[alexgleith]

This is fantastic, I'm going to try to find some time to test it.

🚀

[vincentsarago]

I don't think we want to have statistics and preview implemented because it will be really inefficient for now

[TomAugspurger]

Thanks a ton for working on this @vincentsarago!

[TomAugspurger]

It'd be helpful (for me at least) to list out all the restrictions we have on this DataArray:

1. Can be >2d, but requires that the last two dimensions represent longitude / latitude?
2. Do we care what the names of the lat / lon / x / y dimensions are?
3. Must have a CRS as discovered by .rio.crs
4. ...

Any others we know of?

[vincentsarago]

Do we care what the names of the lat / lon / x / y dimensions are?

reading https://github.com/corteva/rioxarray/blob/dcb26dd0af034dbffd1952bd61d8f62e16b0ee4d/rioxarray/rioxarray.py#L251-L276 it seems that only x,y,longitude,latitude are valid 🤷

[rabernat]

Or CF coordinates with appropriate standard names.

[TomAugspurger]

In case it's helpful for testing, I put a Zarr store up at "https://pangeo.blob.core.windows.net/pangeo-public/daymet-rio-tiler/na-wgs84.zarr/" that's been reprojected to WGS84 and chunked into 256x256 pixel tiles. There's a single time slice, and one variable tmax.

It's Azure's West Europe region.

import xarray
import matplotlib.pyplot as plt
from rio_tiler.io.xarray import XarrayReader

ds = xarray.open_dataset(
    "https://pangeo.blob.core.windows.net/pangeo-public/daymet-rio-tiler/na-wgs84.zarr/",
    engine="zarr",
    decode_coords="all",
    consolidated=True,
)
da = ds["tmax"]

with XarrayReader(da) as dst:
    img = dst.tile(1, 1, 2)

plt.imshow(img.data_as_image());

[geospatial-jeff]

🚀

[vincentsarago]

🙏 thanks everyone

I'm merging this to the main PR and will work on updating the documentation before publishing a new version of rio-tiler with the new XarrayReader.

[TomAugspurger]

Thanks Vincent!

[vietnguyengit]

Hi @TomAugspurger , I had my Conda env setup (Python 3.9.13) with the latest rio-tiler but not sure why I'm unable to import it to the notebook while importing to Python CLI has no issues. Screenshots below:

Successfully imported:

Same Conda env but unable to import to the notebook:

Thanks Tom.

[kylebarron]

This isn't in a released package yet. You can try 4.0.0a0 if you want

[vietnguyengit]

Perfect! Thanks @kylebarron

[vietnguyengit]

Hi @kylebarron @TomAugspurger again,

I wonder why the axes are always 0-250 (my Zarr store is SST ocean data around Australia's region)

and also for tile settings, I wonder if there is somewhere document on how to pick the numbers. I can only get it works with (1,1,1)

When I ran my Zarr store with the XarrayReader fastapi server, the json returned looks like this :

{
"tilejson": "2.1.0",
"name": "xarray",
"version": "1.0.0",
"scheme": "xyz",
"tiles": [
"[http://127.0.0.1:8000/tiles/{z}/{x}/{y}?url=s3%3A%2F%2Fimos-data-pixeldrill%2Fvhnguyen%2Fplayground%2Fmulti-years&variable=sea_surface_temperature](http://127.0.0.1:8000/tiles/%7Bz%7D/%7Bx%7D/%7By%7D?url=s3%3A%2F%2Fimos-data-pixeldrill%2Fvhnguyen%2Fplayground%2Fmulti-years&variable=sea_surface_temperature)"
],
"minzoom": 1,
"maxzoom": 5,
"bounds": [
70.00000213595068,
-69.99999786350668,
190.00000549344387,
19.999999770855315
],
"center": [
130.00000381469727,
-24.999999046325684,
1
]
}

Thank you.

[vincentsarago]

@vietnguyengit is it possible that your X coordinates goes from 0 to 360?

[vietnguyengit]

Thanks @vincentsarago

Actually, I understood what those numbers are after reading the docs here: https://cogeotiff.github.io/rio-tiler/api/rio_tiler/io/xarray/

And the results so far are sweet! 🎉

• Interactive tools helping me pick XY coordinates values: https://www.maptiler.com/google-maps-coordinates-tile-bounds-projection/

Sample results:

with random XY coordinates around Australia's coastal regions

zoom level 2

zoom level 4

zoom level 5

Ping @alexgleith

[alexgleith]

Hi @vietnguyengit, nice work! I'll come find you today, let's see if we can get this showing on a webmap... maybe with a timeslider too!

[yoninachmany]

In case it's helpful for testing, I put a Zarr store up at "https://pangeo.blob.core.windows.net/pangeo-public/daymet-rio-tiler/na-wgs84.zarr/" that's been reprojected to WGS84 and chunked into 256x256 pixel tiles. There's a single time slice, and one variable tmax.

It's Azure's West Europe region.

import xarray
import matplotlib.pyplot as plt
from rio_tiler.io.xarray import XarrayReader

ds = xarray.open_dataset(
    "https://pangeo.blob.core.windows.net/pangeo-public/daymet-rio-tiler/na-wgs84.zarr/",
    engine="zarr",
    decode_coords="all",
    consolidated=True,
)
da = ds["tmax"]

The DataArray has nan values – can they be masked out in tile endpoints or does nodata handling need to be added to XarrayReader.tile?

da
<xarray.DataArray 'tmax' (time: 1, y: 3728, x: 17268)>
array([[[nan, nan, ..., nan, nan],
        [nan, nan, ..., nan, nan],
        ...,
        [nan, nan, ..., nan, nan],
        [nan, nan, ..., nan, nan]]], dtype=float32)

Thanks!

[vincentsarago]

@yoninachmany we didn't implement mask in https://github.com/cogeotiff/rio-tiler/blob/main/rio_tiler/io/xarray.py#L239
I'm not familiar enough with rio-xarray but I feel that it should be possible to create a mask and add it to the ImageData output

Would you be able to open an issue about this ?🙏