Difference between xarray.Dataset.crs and xarray.Dataset.rio.crs #456

geopanda1 · 2022-01-17T11:24:38Z

geopanda1
Jan 17, 2022

Hi,
I have two xarray datasets (Sentinel-2 and MODIS imagery) and I am using rio.reproject_match to reproject MODIS data to the same shape and resolution of the Sentinel-2 data.

Description of the two input datasets:

print(ds_s2) returns

<xarray.Dataset>
Dimensions: (time: 2, y: 1074, x: 862)
Coordinates:

time (time) datetime64[ns] 2017-07-03T10:10:21 2017-07-06T10:20:21

y (y) float64 5.216e+06 5.216e+06 ... 5.205e+06 5.205e+06

x (x) float64 6.531e+05 6.531e+05 ... 6.617e+05 6.617e+05
spatial_ref int32 32632
Data variables:
green (time, y, x) uint16 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
Attributes:
crs: EPSG:32632
grid_mapping: spatial_ref

print(ds_modis) returns

<xarray.Dataset>
Dimensions: (time: 10, y: 25, x: 22)
Coordinates:

time (time) datetime64[ns] 2017-07-01T12:00:00 ... 2017-07-10T12:...

y (y) float64 5.216e+06 5.216e+06 ... 5.206e+06 5.205e+06

x (x) float64 8.365e+05 8.37e+05 ... 8.458e+05 8.462e+05
spatial_ref int32 0
Data variables:
snow_cover (time, y, x) uint8 255 0 0 255 255 255 ... 255 255 255 255 255
Attributes:
crs: PROJCRS["MODIS Sinusoidal",BASEGEOGCRS["WGS 84",DATUM["Wor...
grid_mapping: spatial_ref

print(ds_modis.crs) returns

PROJCRS["MODIS Sinusoidal",BASEGEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["unnamed",METHOD["Sinusoidal"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]

Reprojection

then I run ds_modis_reprojected = ds_modis.rio.reproject_match(ds_s2) and check the crs again.
when I run print(ds_modis_reprojected.crs) I get the following output:

PROJCRS["MODIS Sinusoidal",BASEGEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["unnamed",METHOD["Sinusoidal"],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]

running print(ds_modis_reprojected.rio.crs) returns

EPSG:32632

running print(ds_modis_reprojected) returns

<xarray.Dataset>
Dimensions: (x: 862, y: 1074, time: 10)
Coordinates:

x (x) float64 6.531e+05 6.531e+05 ... 6.617e+05 6.617e+05

y (y) float64 5.216e+06 5.216e+06 ... 5.205e+06 5.205e+06

time (time) datetime64[ns] 2017-07-01T12:00:00 ... 2017-07-10T12:...
spatial_ref int64 0
Data variables:
snow_cover (time, y, x) uint8 255 255 255 255 255 ... 255 255 255 255 0
Attributes:
crs: PROJCRS["MODIS Sinusoidal",BASEGEOGCRS["WGS 84",DATUM["Wor...
grid_mapping: spatial_ref

Question

Obviously, the reprojected dataset seems to be correct in terms of the new dimensions and ds_modis_reprojected.rio.crs represents the correct crs after reprojection. ds_modis_reprojected.crs however still holds the original crs. Is this expected? Is the latter attribute simply not changed by rioxarray?