xr.rio.to_raster() results in gdalbuildvrt does not support positive NS resolution read error in GDAL

[msteckle]

Hello all,

Issue
I think this is a pretty simple question with an easy solution that I'm missing. I am trying to export an xarray dataset using xr.rio.to_raster which is successful, but I have two issues.

1. The outputted .tif cannot be read by gdal when creating a virtual raster (I'm stitching together a LOT of evenly spaced tifs). I get the error gdalbuildvrt does not support positive NS resolution, so I have to do a vaguely annoying workaround to re-project after I run my python script: gdalwarp -overwrite -t_srs EPSG:4326 output.tif corrected_output.tif Then, I can create my virtual raster successfully with the corrected tifs.
2. The CRS information was not saved, which I suppose may be causing issue 1?
   There has to be a way to make sure the pixel size (North-South resolution) is (+,-) instead of (+,+) and to export the CRS information that I'm messing up.

Python Code

# the coords var is a pandas df with EPSG:4326 x and y columns representing lon and lat
# the data var is a 1D array of data values
results = coords.copy()
results['data'] = data
results['band'] = 1 # fill in some value for band

# write 
results_xr = xr.Dataset.from_dataframe(results.set_index(['band', 'y', 'x']))
xr_band = results_xr.isel(band=0) # gets first index band (I only have 1 band)
xr_band.rio.write_crs('EPSG:4326', inplace=True) # I thought this would set my CRS...
xr_band.rio.to_raster("output.tif")

Bash Code

RASTERS=path/to/*.tif
vrt_name=vrtname.vrt
gdalbuildvrt $vrt_name $RASTERS

> gdalbuildvrt does not support positive NS resolution

gdalinfo of output image

(conda env) [user@machine]$ gdalinfo output.tif

> Driver: GTiff/GeoTIFF
> Files: output.tif
> Size is 900, 900
> Origin = (-164.607292662061127,64.355306954322515)
> Pixel Size = (0.000179663056824,0.000179663056824) # <<< The culprit. Also, no CRS info.
> Image Structure Metadata:
>   INTERLEAVE=BAND
> Corner Coordinates:
> Upper Left  (-164.6072927,  64.3553070) 
> Lower Left  (-164.6072927,  64.5170037) 
> Upper Right (-164.4455959,  64.3553070) 
> Lower Right (-164.4455959,  64.5170037) 
> Center      (-164.5264443,  64.4361553) 
> Band 1 Block=900x1 Type=Float64, ColorInterp=Gray
>   Description = datavariable

gdalinfo of re-projected image, which is a step I didn't want to do

(conda env) [user@machine]$ gdalinfo corrected_output.tif

> Driver: GTiff/GeoTIFF
> Files: corrected_output.tif
> Size is 900, 900
> Coordinate System is:
> GEOGCRS["WGS 84", # CRS is now here, woohoo
>     ENSEMBLE["World Geodetic System 1984 ensemble",
>         MEMBER["World Geodetic System 1984 (Transit)"],
>         MEMBER["World Geodetic System 1984 (G730)"],
>         MEMBER["World Geodetic System 1984 (G873)"],
>         MEMBER["World Geodetic System 1984 (G1150)"],
>         MEMBER["World Geodetic System 1984 (G1674)"],
>         MEMBER["World Geodetic System 1984 (G1762)"],
>         MEMBER["World Geodetic System 1984 (G2139)"],
>         ELLIPSOID["WGS 84",6378137,298.257223563,
>             LENGTHUNIT["metre",1]],
>         ENSEMBLEACCURACY[2.0]],
>     PRIMEM["Greenwich",0,
>         ANGLEUNIT["degree",0.0174532925199433]],
>     CS[ellipsoidal,2],
>         AXIS["geodetic latitude (Lat)",north,
>             ORDER[1],
>             ANGLEUNIT["degree",0.0174532925199433]],
>         AXIS["geodetic longitude (Lon)",east,
>             ORDER[2],
>             ANGLEUNIT["degree",0.0174532925199433]],
>     USAGE[
>         SCOPE["Horizontal component of 3D system."],
>         AREA["World."],
>         BBOX[-90,-180,90,180]],
>     ID["EPSG",4326]]
> Data axis to CRS axis mapping: 2,1
> Origin = (-165.577473168910188,64.517003705464035)
> Pixel Size = (0.000179663056824,-0.000179663056824) # Fixed
> Metadata:
>   AREA_OR_POINT=Area
> Image Structure Metadata:
>   INTERLEAVE=BAND
> Corner Coordinates:
> Upper Left  (-165.5774732,  64.5170037) (165d34'38.90"W, 64d31' 1.21"N)
> Lower Left  (-165.5774732,  64.3553070) (165d34'38.90"W, 64d21'19.11"N)
> Upper Right (-165.4157764,  64.5170037) (165d24'56.80"W, 64d31' 1.21"N)
> Lower Right (-165.4157764,  64.3553070) (165d24'56.80"W, 64d21'19.11"N)
> Center      (-165.4966248,  64.4361553) (165d29'47.85"W, 64d26'10.16"N)
> Band 1 Block=900x1 Type=Float64, ColorInterp=Gray
>   Description = fcover

---

And here is some system information you might want to know. I'm SSH'ed into a remote system while doing this.
Python Environment

sys.version
> '3.11.5 | packaged by conda-forge | (main, Aug 27 2023, 03:34:09) [GCC 12.3.0]'

import rioxarray as rxr
rxr.__version__
> '0.15.0'

import xarray as xr
xr.__version__
> '2023.10.1'

GDAL and Conda

gdalinfo --version
> GDAL 3.7.1, released 2023/07/06

conda --version
> conda 4.9.2

Thanks!
Morgan

[snowman2]

does not support positive NS resolution

See: pydata/xarray#6695

xr_band = xr_band.isel({xr_band.rio.y_dim: slice(None, None, -1)})

The CRS information was not saved

Not sure about this one. Would need a re-producible example. https://stackoverflow.com/help/minimal-reproducible-example

[RichardScottOZ]

Possibly related? #692