Reproject indexing errors #762
Replies: 4 comments 7 replies
-
|
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:
|
Beta Was this translation helpful? Give feedback.
-
|
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 |
Beta Was this translation helpful? Give feedback.
-
|
For the 2.2km UKCP dataset on a rotated pole grid, If you force rioxarray to use As you can see below: #762 (comment) |
Beta Was this translation helpful? Give feedback.
-
|
@snowman2 looks like this data is using location arrays for geo referencing. Is this something Would be nice to agree on a format for dealing with location arrays. |
Beta Was this translation helpful? Give feedback.
-
Beta Was this translation helpful? Give feedback.
-
|
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 |
Beta Was this translation helpful? Give feedback.
-
|
So: I think I've gotten a lot further by removing the 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: nanInitially I thought that was a success (and note the Then my colleague @andrewphilipsmith checked that via Hoping the following screenshots make this clear: 
|
Beta Was this translation helpful? Give feedback.
-
|
@griff-rees Did you get anywhere with this? I have the exact same issue 😅 Same files too! |
Beta Was this translation helpful? Give feedback.
-
|
Pls see my comment above #762 (reply in thread) - the 1D |
Beta Was this translation helpful? Give feedback.
-
|
Best bet if you want to use Rasterio is probably to try the new Else, you can use a regridding software like xESMF, as mentioned here #209 (comment): 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. |
Beta Was this translation helpful? Give feedback.
-
|
Btw, this |
Beta Was this translation helpful? Give feedback.
Uh oh!
There was an error while loading. Please reload this page.
Uh oh!
There was an error while loading. Please reload this page.
-
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:My guess is I need to figure the exact coordinate structure to solve the
undefinedcomponents, and I wonder if it relates toEURO-CORDEXand/orcartopyRotatedGeodetic. While I try to figure that (at which point I assume I could then usedataset.rio.write_crs(), this is how far I've got via other routes:I've also tried using
set_spatial_dimsbut that eventually leads to the same1-dimensionalerror.Update
I think this relates to #772. I've ended up going between
rioxarrayandosgeo.gdalacross various steps, hope to post below if I reach a solution.Beta Was this translation helpful? Give feedback.
All reactions