I have a couple GeoTiffs where the "no data" value is set to 32767. I want to set those values to NaN and then use nearest neighbor interpolation to resample the data and match the projection and resolution of the other GeoTiff. To do this, I use rio.reproject_match. The problem is that when I reproject the layer, the sum of all items returns 0 which I don't think is correct. Can someone find what I'm doing wrong here?
import numpy as np
import rioxarray
import xarray
from rasterio.enums import Resampling
# Source Raster
xds = rioxarray.open_rasterio(LS)
# Set no data value 32767 to NaN
xds = xds.where(xds != 32767, np.nan)
xds.rio.write_nodata(np.nan, inplace=True)
print(xds)
# <xarray.DataArray (band: 1, y: 1006, x: 950)>
# 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]]])
# Coordinates:
# * band (band) int64 1
# * x (x) float64 -4.56e+03 -4.53e+03 ... 2.388e+04 2.391e+04
# * y (y) float64 2.029e+06 2.029e+06 ... 1.999e+06 1.998e+06
# spatial_ref int64 0
# Attributes: (12/13)
# AREA_OR_POINT: Area
# DataType: Generic
# RepresentationType: ATHEMATIC
# STATISTICS_COVARIANCES: 20807216.58662666
# STATISTICS_MAXIMUM: 11111
# STATISTICS_MEAN: 2561.3518011412
# ... ...
# STATISTICS_SKIPFACTORX: 1
# STATISTICS_SKIPFACTORY: 1
# STATISTICS_STDDEV: 4561.4928024307
# _FillValue: nan
# scale_factor: 1.0
# add_offset: 0.0
print_raster(xds)
# Source Raster:
# ----------------
# shape: (1006, 950)
# resolution: (30.0, -30.0)
# bounds: (-4575.0, 1998465.0, 23925.0, 2028645.0)
# sum: 2284131573.0
# CRS: EPSG:9822
# Target Raster
xds_match = rioxarray.open_rasterio(precip)
xds_match = xds_match.where(xds_match != 32767, np.nan)
xds_match.rio.write_nodata(np.nan, inplace=True)
print_raster(xds_match)
# Raster to Match:
# ----------------
# shape: (2898, 2772)
# resolution: (10.0, -10.0)
# bounds: (747720.0, 4543380.0, 775440.0, 4572360.0)
# sum: 7517694.093363524
# CRS: EPSG:32614
# Resample and reproject the source raster to match the target raster
xds_repr_match = xds.rio.reproject_match(xds_match, resampling=Resampling.nearest, nodata=np.nan)
print_raster(xds_repr_match)
# Reprojected Raster:
# -------------------
# shape: (2898, 2772)
# resolution: (10.0, -10.0)
# bounds: (747720.0, 4543380.0, 775440.0, 4572360.0)
# sum: 0.0
# CRS: EPSG:32614
I don't expect the sum to match the target raster sum, but I do expect it to be more than 0.
I've read this thread, but they're talking about reproject instead of reproject_match mostly. How to deal with Nan value when using rioxarray rio.reproject()?. Nonetheless, I built a little test case using reproject with some dummy data just to see if I understood how sum works in conjunction with NaN values and it does it properly (doesn't return 0):
import xarray as xr
import rioxarray as rxr
import numpy as np
from rasterio.enums import Resampling
data = [[1, 2, 3, 4], [5, 32767, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]]
x = np.arange(100, 116, 4)
y = np.arange(36, 20, -4)
foo = xr.DataArray(data, coords=[y, x], dims=['y', 'x'])
foo = foo.rio.write_crs('EPSG:4326')
foo = foo.where(foo != 32767, np.nan)
foo.rio.write_nodata(np.nan, inplace=True)
foo_resampled = foo.rio.reproject(
foo.rio.crs,
resolution=8,
resampling=Resampling.nearest,
nodata=np.nan
)
print(foo.sum().item())
# 130.0
print(foo_resampled.sum().item())
# 38.0
Why doesn't reproject_match work in the way I expect with np.nan values and how can I fix it?
