corteva / rioxarray

geospatial xarray extension powered by rasterio
https://corteva.github.io/rioxarray
Other
504 stars 80 forks source link

Unable to generate consistent resolution #690

Closed dluks closed 11 months ago

dluks commented 11 months ago

I'm trying to create an ideal synthetic placeholder dataset to use as a target for reprojecting some real datasets using reproject_match. I'm not sure if this is a bug or simply user error/misunderstanding, but it's giving me a minor headache. As much as I try to eliminate float-precision errors with rounding and an even division of the desired resolution into the coordinates, I can't seem to escape them when inspecting them using dataset.rio.resolution().

import numpy as np
import rioxarray as riox
import xarray as xr

# Make an empty dataset with the new extent and resolution
# Define the global extent
xmin, ymin, xmax, ymax = (-180.0, -90.0, 180.0, 90.0)

# Calculate the number of rows and columns based on the desired resolution
res = 0.01
nrows = int((ymax - ymin) / res)
ncols = int((xmax - xmin) / res)

# Generate lat (y) and lon (x) grid and round values to two decimal places
y = np.round(np.linspace(ymax, ymin, nrows), 2)
x = np.round(np.linspace(xmin, xmax, ncols), 2)

# Create the empty DataArray with the desired dimensions and fill with ones
target_grid = xr.DataArray(
    data=np.ones((nrows, ncols)),
    dims=("y", "x"),
    coords={
        "y": y,
        "x": x,
    },
)
target_grid.rio.write_crs("EPSG:4326", inplace=True)

Problem description

Despite doing everything in my power to create a clean grid uncontaminated with float-precision errors, the resulting resolution still comes out imperfect (I would have expected a simple (0.01, -0.01)):

target_grid.rio.resolution()  # (0.010000277785494041, -0.010000555586421468)

Environment Information

rioxarray (0.14.0) deps:
  rasterio: 1.3.7
    xarray: 2023.7.0
      GDAL: 3.7.0
      GEOS: 3.11.2
      PROJ: 9.2.1
 PROJ DATA: /home/lusk/mambaforge/envs/traits/share/proj
 GDAL DATA: /home/lusk/mambaforge/envs/traits/share/gdal

Other python deps:
     scipy: 1.10.1
    pyproj: 3.6.0

System:
    python: 3.9.16 | packaged by conda-forge | (main, Feb  1 2023, 21:39:03)  [GCC 11.3.0]
executable: /home/lusk/mambaforge/envs/traits/bin/python
   machine: Linux-5.13.0-27-generic-x86_64-with-glibc2.31
snowman2 commented 11 months ago

When generating coordinates, it is important to remember that the coordinates represent the center of the grid cell. That means you have to adjust by resolution/2 from the outer bounds.

dluks commented 11 months ago

🤦‍♂️ Of course. I had a feeling it was user error. That solved the issue. Thanks!