Pixel differences between GDAL and Rioaxarray rasters

[nlamkey]

I'm working on creating a raster stack and did a compare and contrast between GDAL and Rioxarray mostly for processing speed and noticed some differences in outputs between to the two.
All I'm doing is opening, reprojecting, resampling and clipping the rasters. Here's my code for Rioxarray:

import os

import numpy as np
import rioxarray as rxr
import xarray as xr
import rasterio as rio
from rasterio.plot import plotting_extent
from rasterio.mask import mask
from rasterio.enums import Resampling

import pyproj
from shapely.geometry import mapping
from osgeo import gdal
import fiona
pyproj.set_use_global_context(True)

prism_path = os.path.join(data_path,
                         "PRISM",
                         "prism_precip",
                          "PRISM_ppt_stable_4kmM3_2017_all_bil",
                          "sum",
                          "PRISM_ppt_stable_4kmM3_2017_bil.bil")

# Open Site Boundary
site_bound = gpd.read_file(site_boundary_path)
site_bound.crs
with fiona.open(site_boundary_path, "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]

# Use rioxarray to open, resample, reproject and clip raster data
prism = rxr.open_rasterio(prism_path,
                                 masked=True)
prism_rep = prism.rio.reproject(site_bound.crs,
                                shape=(6797,12274),
                                resampling=Resampling.nearest)

prism_clip = prism_rep.rio.clip(geometries=shapes,
                            all_touched=True,
                            from_disk=True).squeeze()

# Plot the data
precip_extents = plotting_extent(prism_clip.values,
                                prism_clip.rio.transform())

f, ax = plt.subplots()

im = ax.imshow(prism_clip,
               extent=precip_extents)
# site_bound.plot(ax=ax)
plt.show()

Which returns this plot:

Here's my code for GDAL:

prism_precip = gdal.Open(prism_path)
array = prism_precip.GetRasterBand(1).ReadAsArray()

# Reproject, resample and clip the data
prism_precip_rep = gdal.Warp('prism_pre_rep.tif', prism_precip, dstSRS=site_bound.crs)

prism_precip_resample = gdal.Warp('prism_pre_res.tif',prism_precip_rep, xRes=500, yRes =500,
                                  resampleAlg = "near")
                             
prism_precip_clip = gdal.Warp('prism_clip.tif', prism_precip_resample,cropToCutline=True, cutlineDSName=site_boundary_path,
                              dstAlpha=True)

# Export the data
band = prism_precip_clip.GetRasterBand(1)
gt = prism_precip_clip.GetGeoTransform()
proj = prism_precip_clip.GetProjection()
array = band.ReadAsArray()

driver = gdal.GetDriverByName("GTiff")
driver.Register()
outds = driver.Create("prism_precip_clip.tif",xsize=array.shape[1],ysize=array.shape[0],
                      bands=1,eType=gdal.GDT_Int32)

outds.SetGeoTransform(gt)
outds.SetProjection(proj)
outband = outds.GetRasterBand(1)
outband.WriteArray(array)
outband.FlushCache()

outband = None
outds= None

# Plot the data
prism = gdal.Open('prism_clip.tif')
band = prism.GetRasterBand(1)
array = band.ReadAsArray()
f, ax = plt.subplots()
ax.imshow(array)

Which returns this plot:

As you can see the pixels in the rioxarray plot look tilted t while the GDAL pixels look square. I wanted to get an opinion on why this is the case and if it really matters for a raster stack or if there's something I'm doing in my code that's causing this.

[snowman2]

@nlamkey,

1. The code you have is not re-producible as you have a lot of missing pieces.
2. It appears the code for rioxarray and GDAL are not doing the exact same thing (some differneces there). Would you mind updating the code snippets so they are doing the same thing?

Things that are sticking out to me that may be helpful:

1. Why are you re-projecting before clipping?
2. from_disk=True only works when performing the clip right after opening the file. It does not work if you re-project first. ref

[nlamkey]

Thanks for the reply. I guess I'm not aware that they're not doing the same thing? Can you point out where they are different?

I am reprojecting before clipping because the clip function requires a site_bound.crs argument that puts the raster and vector layer into the same crs while the clip takes place and I couldn't find a good way to resample without using reproject(). So I didn't want to have to reproject the data twice.

Also doing the clip first and the reproject second gave me different cell sizes than doing it the other way around.

[snowman2]

With rioxarray, you pass in shape=(6797,12274), and with GDAL you pass in xRes=500, yRes =500,. rio.reproject has a resolution argument you can pass in to make them the same.

[snowman2]

Also: https://corteva.github.io/rioxarray/stable/examples/reproject.html#Reproject-Large-Rasters-with-Virtual-Warping

with fiona.open(site_boundary_path, "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]
    crs = shapefile.crs

# Use rioxarray to open, resample, reproject and clip raster data
with rasterio.open(prism_path) as src, rasterio.vrt.WarpedVRT(src, crs=crs) as vrt:
    prism = rxr.open_rasterio(vrt, masked=True).rio.clip(
        geometries=shapes,
        all_touched=True,
        from_disk=True,
    ).squeeze()