Reproject indexing errors #762

griff-rees · 2024-03-27T16:02:15Z

griff-rees
Mar 27, 2024

Thanks for the library. I'm attempting to process https://catalogue.ceda.ac.uk/uuid/d5822183143c4011a2bb304ee7c0baf7, reprojected to EPSG:27700. The (I think) related examples I've found thus far include:

>>> from pathlib import Path
>>> from xarray import open_dataset, Dataset
>>> import rioxarray

>>> url: Path = Path('tests/data/tasmax_rcp85_land-cpm_uk_2.2km_01_day_19801201-19811130.nc')

>>> dataset: Dataset = open_dataset(url, decode_coords="all")
>>> dataset
<xarray.Dataset> Size: 427MB
Dimensions:                     (ensemble_member: 1, time: 360,
                                 grid_latitude: 606, grid_longitude: 484,
                                 bnds: 2)
Coordinates: (12/14)
    rotated_latitude_longitude  int32 4B -2147483647
  * ensemble_member             (ensemble_member) int32 4B 1
  * time                        (time) object 3kB 1980-12-01 12:00:00 ... 198...
    time_bnds                   (time, bnds) object 6kB ...
  * grid_latitude               (grid_latitude) float64 5kB -4.683 ... 8.063
    grid_latitude_bnds          (grid_latitude, bnds) float64 10kB ...
    ...                          ...
    ensemble_member_id          (ensemble_member) |S27 27B ...
    latitude                    (grid_latitude, grid_longitude) float64 2MB 4...
    longitude                   (grid_latitude, grid_longitude) float64 2MB ...
    month_number                (time) int32 1kB ...
    year                        (time) int32 1kB 1980 1980 1980 ... 1981 1981
    yyyymmdd                    (time) |S64 23kB b'19801201                  ...
Dimensions without coordinates: bnds
Data variables:
    tasmax                      (ensemble_member, time, grid_latitude, grid_longitude) float32 422MB ...
Attributes: (12/15)
    collection:      land-cpm
    contact:         ukcpproject@metoffice.gov.uk
    creation_date:   2021-05-11T14:06:30
    domain:          uk
    frequency:       day
    institution:     Met Office Hadley Centre (MOHC), FitzRoy Road, Exeter, D...
    ...              ...
    resolution:      2.2km
    scenario:        rcp85
    source:          UKCP18 realisation from a set of 12 convection-permittin...
    title:           UKCP18 land projections - 2.2km convection-permitting cl...
    version:         v20210615
    Conventions:     CF-1.7
>>> dataset.rio.crs
CRS.from_wkt('GEOGCRS["undefined",BASEGEOGCRS["undefined",DATUM["undefined",ELLIPSOID["undefined",6371229,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["undefined",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],DERIVINGCONVERSION["Pole rotation (netCDF CF convention)",METHOD["Pole rotation (netCDF CF convention)"],PARAMETER["Grid north pole latitude (netCDF CF convention)",37.5,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["Grid north pole longitude (netCDF CF convention)",177.5,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["North pole grid longitude (netCDF CF convention)",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]')

My guess is I need to figure the exact coordinate structure to solve the undefined components, and I wonder if it relates to EURO-CORDEX and/or cartopy RotatedGeodetic. While I try to figure that (at which point I assume I could then use dataset.rio.write_crs(), this is how far I've got via other routes:

>>> from typing import Final
>>> FINAL_CRS: Final[str] = 'EPSG:27700'
>>> dataset.rio.reproject(FINAL_CRS)
*** rioxarray.exceptions.TooManyDimensions: Only 2D and 3D data arrays supported. Data variable: tasmax
>>> xr_time_series.squeeze().rio.reproject(UK_SPATIAL_PROJECTION)
*** ValueError: IndexVariable objects must be 1-dimensional

I've also tried using set_spatial_dims but that eventually leads to the same 1-dimensional error.

Update

I think this relates to #772. I've ended up going between rioxarray and osgeo.gdal across various steps, hope to post below if I reach a solution.

snowman2 · 2024-03-27T19:55:30Z

snowman2
Mar 27, 2024
Maintainer

If I had a file, this is what I would try:

dataset[["rotated_latitude_longitude", "grid_latitude", "grid_longitude", "tasmax"]].isel(
    ensemble_member=0,
    time=0,
).rio.set_spatial_dims(
    x_dim="grid_longitude",
    y_dim="grid_latitude",
    inplace=True,
).rio.reproject("'EPSG:27700")

Steps:

Only selected the variables with the information needed
Convert grid from 4D to 2D by selecting a specific time slice

2 replies

griff-rees Mar 28, 2024
Author

Thanks @snowman2! Sadly same error but will keep trying variations:

>>> dataset[["rotated_latitude_longitude", "grid_latitude", "grid_longitude", "tasmax"]].isel(
...        ensemble_member=0, 
...        time=0,
...    ).rio.set_spatial_dims(
...        x_dim="grid_longitude",
...        y_dim="grid_latitude",
...        inplace=True,
...    ).rio.reproject("EPSG:27700")
*** ValueError: IndexVariable objects must be 1-dimensional
>>> dataset[["rotated_latitude_longitude", "grid_latitude", "grid_longitude", "tasmax", "bnds"]].isel(
...        ensemble_member=0, 
...        time=0, 
...        bnds=0,
...   ).rio.set_spatial_dims(
...        x_dim="grid_longitude", 
...        y_dim="grid_latitude", 
...        inplace=True,
...   ).rio.reproject("EPSG:27700")
*** ValueError: IndexVariable objects must be 1-dimensional

j08lue Sep 5, 2024

For the 2.2km UKCP dataset on a rotated pole grid, rio.set_spatial_dims(x_dim="grid_longitude", y_dim="grid_latitude") unfortunately is wrong, because those 1D (rectilinear) coordinate axes are incorrect wrt the actual (curvilinear) projection. Only the 2D latitude and longitude arrays are correct to feed into reprojection - or a specification of the rotated grid.

If you force rioxarray to use grid_longitude and grid_latitude, you get a result, but it will be wrongly georeferenced, especially visible at the edges.

As you can see below: #762 (comment)

Kirill888 · 2024-03-29T00:58:31Z

Kirill888
Mar 29, 2024

@snowman2 looks like this data is using location arrays for geo referencing. Is this something rioxarray can represent? In odc-geo I had to work around by creating a GCPGeoBox from samples read from location array when dealing with data like that.

Would be nice to agree on a format for dealing with location arrays.

1 reply

snowman2 Mar 29, 2024
Maintainer

See: rasterio/rasterio#2884

griff-rees · 2024-05-24T22:05:44Z

griff-rees
May 24, 2024
Author

Thanks again @Kirill888 and @snowman2. FYI: looking another route: https://gis.stackexchange.com/questions/312084/transform-coordinate-from-rotated-lat-lon-to-normal-lat-lon-wgs84

0 replies

griff-rees · 2024-06-24T11:35:51Z

griff-rees
Jun 24, 2024
Author

So: I think I've gotten a lot further by removing the ensemble_member dim via the .sel() method, but run into an issue with the projection

The current process:

>>> from pathlib import Path
>>> from xarray import open_dataset, Dataset, DataArray
>>> import numpy as np
>>> import rioxarray

>>> cpm_path: Path = Path('tests/data/local-cache/tasmax_rcp85_land-cpm_uk_2.2km_01_day_19801201-19811130.nc')
>>> dataset: Dataset = open_dataset(cpm_path, decode_coords="all")
>>> dataset
<xarray.Dataset> Size: 427MB
Dimensions:                     (ensemble_member: 1, time: 360,
                                 grid_latitude: 606, grid_longitude: 484,
                                 bnds: 2)
Coordinates: (12/14)
    rotated_latitude_longitude  int32 4B -2147483647
  * ensemble_member             (ensemble_member) int32 4B 1
  * time                        (time) object 3kB 1980-12-01 12:00:00 ... 198...
    time_bnds                   (time, bnds) object 6kB ...
  * grid_latitude               (grid_latitude) float64 5kB -4.683 ... 8.063
    grid_latitude_bnds          (grid_latitude, bnds) float64 10kB ...
    ...                          ...
    ensemble_member_id          (ensemble_member) |S27 27B ...
    latitude                    (grid_latitude, grid_longitude) float64 2MB 4...
    longitude                   (grid_latitude, grid_longitude) float64 2MB ...
    month_number                (time) int32 1kB ...
    year                        (time) int32 1kB 1980 1980 1980 ... 1981 1981
    yyyymmdd                    (time) |S64 23kB b'19801201                  ...
Dimensions without coordinates: bnds
Data variables:
    tasmax                      (ensemble_member, time, grid_latitude, grid_longitude) float32 422MB ...
Attributes: (12/15)
    collection:      land-cpm
    contact:         ukcpproject@metoffice.gov.uk
    creation_date:   2021-05-11T14:06:30
    domain:          uk
    frequency:       day
    institution:     Met Office Hadley Centre (MOHC), FitzRoy Road, Exeter, D...
    ...              ...
    resolution:      2.2km
    scenario:        rcp85
    source:          UKCP18 realisation from a set of 12 convection-permittin...
    title:           UKCP18 land projections - 2.2km convection-permitting cl...
    version:         v20210615
    Conventions:     CF-1.7
>>> dataset.rio.crs
CRS.from_wkt('GEOGCRS["undefined",BASEGEOGCRS["undefined",DATUM["undefined",ELLIPSOID["undefined",6371229,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["undefined",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],DERIVINGCONVERSION["Pole rotation (netCDF CF convention)",METHOD["Pole rotation (netCDF CF convention)"],PARAMETER["Grid north pole latitude (netCDF CF convention)",37.5,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["Grid north pole longitude (netCDF CF convention)",177.5,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["North pole grid longitude (netCDF CF convention)",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]')
>>> without_ensemble: DataArray = dataset.sel(ensemble_member=1)
>>> coords: dict[str, DataArray] = {
...   'time': dataset['time'], 
...   'grid_latitude': dataset['grid_latitude'], 
...   'grid_longitude': dataset['grid_longitude']
... }
>>> for_reprojection = DataArray(data=without_ensemble.to_numpy(), coords=coords, name="tasmax")
>>> for_reprojection.rio.write_crs(input_crs=dataset.rio.crs, inplace=True)
>>> for_reprojection.rio.reproject(dst_crs="EPSG:27700", nodata=np.nan)
<xarray.DataArray 'tasmax' (time: 360, y: 653, x: 529)> Size: 497MB
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
...
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
Coordinates:
  * x            (x) float64 4kB -3.129e+05 -3.107e+05 ... 8.465e+05 8.487e+05
  * y            (y) float64 5kB 1.197e+06 1.195e+06 ... -2.353e+05 -2.375e+05
  * time         (time) object 3kB 1980-12-01 12:00:00 ... 1981-11-30 12:00:00
    spatial_ref  int64 8B 0
Attributes:
    _FillValue:  nan

Initially I thought that was a success (and note the nan values are expected on the edges).

Then my colleague @andrewphilipsmith checked that via QGIS and found alignment issues. Any suggestions appreciated.

Hoping the following screenshots make this clear:

A map of the UK with red circles around coast points where a blue border doesn't align with rainall measurements

A map of a northern coast of Scotland with a red layer out of alignment with a blue border line

4 replies

import xarray as xr
import xesmf as xe
import numpy as np
import rioxarray

ds = xr.open_dataset("curvilinear.nc")

dlon = dlat = 0.01
target_grid = {
    "lon": np.arange(ds["longitude"].min(), ds["longitude"].max() + dlon, dlon, dtype="float32"), 
    "lat": np.arange(ds["latitude"].min(), ds["latitude"].max() + dlat, dlat, "float32")
}

regridder = xe.Regridder(ds, target_grid, "bilinear", unmapped_to_nan=True)
ds_regridded = regridder(ds)

da_regridded = ds_regridded[VARIABLE_NAME]

da_regridded.rio.write_crs("epsg:4326", inplace=True)
da_regridded.rio.set_spatial_dims("lon", "lat", inplace=True)
da_regridded.rio.to_raster("rectilinear.tif")

The above works for the UKCP 2.2km file and produces a correctly geo-referenced GeoTIFF.

j08lue Sep 5, 2024

Btw, this src_geoloc_array GDAL method that @snowman2 got into Rasterio (coming in Rasterio 1.4) also works for the UKCP 2.2km UK dataset: #209 (comment).

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Reproject indexing errors #762

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 4 comments 7 replies

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Select a reply

Reproject indexing errors #762

griff-rees Mar 27, 2024

Update

Replies: 4 comments · 7 replies

snowman2 Mar 27, 2024 Maintainer

griff-rees Mar 28, 2024 Author

j08lue Sep 5, 2024

Kirill888 Mar 29, 2024

snowman2 Mar 29, 2024 Maintainer

griff-rees May 24, 2024 Author

griff-rees Jun 24, 2024 Author

ciaransweet Sep 4, 2024

j08lue Sep 5, 2024

j08lue Sep 5, 2024

j08lue Sep 5, 2024

griff-rees
Mar 27, 2024

Replies: 4 comments 7 replies

snowman2
Mar 27, 2024
Maintainer

griff-rees Mar 28, 2024
Author

Kirill888
Mar 29, 2024

snowman2 Mar 29, 2024
Maintainer

griff-rees
May 24, 2024
Author

griff-rees
Jun 24, 2024
Author