Parsing coordinates

[norlandrhagen]

Hey there!

I'm trying to combine some tiffs with open_virtual_mfdataset, but when reading these ESA tiffs, they don't have any coordinate information or CRS. Is there any functionality for inferring CRS from the tiff?

import obstore
from virtualizarr import open_virtual_dataset
from virtualizarr.registry import ObjectStoreRegistry
from virtual_tiff import VirtualTIFF

bucket_url = "s3://esa-worldcover/"
file_url = f"{bucket_url}v100/2020/map/ESA_WorldCover_10m_2020_v100_N24W075_Map.tif"

s3_store = obstore.store.from_url(bucket_url, region="eu-central-1", skip_signature=True)
registry = ObjectStoreRegistry({bucket_url: s3_store})

ds = open_virtual_dataset(
    url=file_url,
    registry=registry,
    parser=VirtualTIFF(ifd=0) 
)
ds

<xarray.Dataset> Size: 1GB
Dimensions:  (y: 36000, x: 36000)
Dimensions without coordinates: y, x
Data variables:
    0        (y, x) uint8 1GB ManifestArray<shape=(36000, 36000), dtype=uint8...

vs rioxarray:

import rioxarray

ds = xr.open_dataset(file_url,engine='rasterio')
print(ds)

<xarray.Dataset> Size: 5GB
Dimensions:      (band: 1, x: 36000, y: 36000)
Coordinates:
  * band         (band) int64 8B 1
  * x            (x) float64 288kB -75.0 -75.0 -75.0 -75.0 ... -72.0 -72.0 -72.0
  * y            (y) float64 288kB 27.0 27.0 27.0 27.0 ... 24.0 24.0 24.0 24.0
    spatial_ref  int64 8B ...
Data variables:
    band_data    (band, y, x) float32 5GB ...

[maxrjones]

thanks again for trying this out @norlandrhagen!

We definitely should provide some convenience functions for adding geospatial coordinates. I sorta side-tracked to try to figure out how to store the required metadata as Zarr attributes (xref https://github.com/zarr-conventions/geo-proj), but now aim to get this library working better. For immediate use, here's how to do the assignment manually using xproj and rasterix. These are all very bleeding edge, so thank you for trying it out and please share again here or in those repos if you notice anything.

from obstore.store import HTTPStore
from virtualizarr import open_virtual_dataset
from virtualizarr.registry import ObjectStoreRegistry
from virtual_tiff import VirtualTIFF
from urllib.parse import urlparse

import xarray as xr
import rioxarray
import zarr
import xproj
import rasterix
from affine import Affine
%matplotlib inline

zarr.config.set({'async.concurrency': 128})

url = "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/18/S/VH/2018/1/S2B_T18SVH_20180116T160106_L2A/TCI.tif"


parsed = urlparse(url)
# Setup and open dataset
store = HTTPStore.from_url(f"{parsed.scheme}://{parsed.netloc}")
registry = ObjectStoreRegistry({url: store})

parser=VirtualTIFF(ifd=0)

vds = open_virtual_dataset(
    url=url,
    registry=registry,
    parser=parser,
)

ms = parser(url=url, registry=registry)
ds = xr.open_zarr(ms, consolidated=False, zarr_format=3)
da_vt = ds["0"].load()
da_vt = da_vt.proj.assign_crs(spatial_ref=f"EPSG:{da_vt.attrs['projected_type']}")
model_pixel_scale = da_vt.attrs['model_pixel_scale']
model_tiepoint = da_vt.attrs['model_tiepoint']
transform = Affine(
    model_pixel_scale[0],  # pixel width
    0.0,                    # rotation
    model_tiepoint[3],      # origin x
    0.0,                    # rotation
    -model_pixel_scale[1],  # pixel height (negative)
    model_tiepoint[4]       # origin y
)
index = rasterix.RasterIndex.from_transform(transform, width=da_vt.sizes['x'], height=da_vt.sizes['y'])
coords = xr.Coordinates.from_xindex(index)
da_vt = da_vt.assign_coords(coords)
da_vt.plot.imshow(rgb="band") 

Here's an example with your data (but I didn't plot due to bad loading performance. I'll look into that)

from virtualizarr import open_virtual_dataset
from virtualizarr.registry import ObjectStoreRegistry
from virtual_tiff import VirtualTIFF
from urllib.parse import urlparse

import obstore
import xarray as xr
import rioxarray
import zarr
import xproj
import rasterix
from affine import Affine
%matplotlib inline

zarr.config.set({'async.concurrency': 128})

bucket_url = "s3://esa-worldcover/"
url = f"{bucket_url}v100/2020/map/ESA_WorldCover_10m_2020_v100_N24W075_Map.tif"

s3_store = obstore.store.from_url(bucket_url, region="eu-central-1", skip_signature=True)
registry = ObjectStoreRegistry({bucket_url: s3_store})

parser=VirtualTIFF(ifd=0)
ms = parser(url=url, registry=registry)
native_data = xr.open_dataset(ms, engine="zarr", consolidated=False, zarr_format=3)["0"]

native_data = native_data.proj.assign_crs(spatial_ref=f"EPSG:{native_data.attrs['geographic_type']}")
model_pixel_scale = native_data.attrs['model_pixel_scale']
model_tiepoint = native_data.attrs['model_tiepoint']
transform = Affine(
    model_pixel_scale[0],  # pixel width
    0.0,                    # rotation
    model_tiepoint[3],      # origin x
    0.0,                    # rotation
    -model_pixel_scale[1],  # pixel height (negative)
    model_tiepoint[4]       # origin y
)
index = rasterix.RasterIndex.from_transform(transform, width=native_data.sizes['x'], height=native_data.sizes['y'])
coords = xr.Coordinates.from_xindex(index)
native_data = native_data.assign_coords(coords)

[norlandrhagen]

Thanks @maxrjones! This is super helpful. Kind of tangential, but have you tried out combining virtual tiffs with raster indexes into a virtual dataset? It seems like the rasterix combining options are a bit more limited.

[maxrjones]

Thanks @maxrjones! This is super helpful. Kind of tangential, but have you tried out combining virtual tiffs with raster indexes into a virtual dataset? It seems like the rasterix combining options are a bit more limited.

Not yet, but I hope to soon. @sharkinsspatial has been working with concatenation via rasterix this week. I think he also ran into some issues

[sharkinsspatial]

@norlandrhagen I'm only working with individual virtual tiffs rather than combining them, so that I can do windowed reads for tiling. I'm hitting a bit of an issue with the how rasterix infers formula based coordinates that I think @dcherian is looking into.

[dcherian]

combining virtual tiffs with raster indexes into a virtual dataset? It seems like the rasterix combining options are a bit more limited.

Do you have an example?