Raster resample to grid with nearest method and nodata vals

[javgs-bd]

Hello. I understand that this line in the current resample_to_grid code: if extrapolate_edges and method != "nearest": results in that all other resampling methods will fill the values for cells that touch raster voids (using the nearest method). The only method that will yield nodatavals for those grid cells is, actually, the nearest method. This means those cells need to be dealt with after the resampling.

I wonder if it would not be more intuitive that the extrapolate_edges option could also be used with the nearest method. Maybe I'm missing something in the resample function usage.

[javgs-bd]

Just a note, I have modified my raster.py file to leave only if extrapolate_edges: in line 443 and it allows me to resample the rasters with the nearest method in my model cells falling outside of the raster extension.

[javgs-bd]

I wonder if the current behavior of extrapolate_edges is the intended one, or if I should open an issue...

[jdhughes-usgs]

@jlarsen-usgs you should probably weigh in on this discussion item

[jlarsen-usgs]

@javgs-bd

There is in fact a small bug in the raster resampling, however it is unrelated to the extrapolate_edges flag. Nearest neighbor, by design, will set the nearest raster value to a grid cell even when the grid cell is outside of the raster boundaries. This falls apart however, when there are no data values (nan values) around an edge/edges of the raster. In the extrapolate edges code we filter out those nan values when they exist, however in the current version of the code when estrapolate_edges=False we were not. I'll be pushing an update today to filter nan values from the raster prior to resampling.

[javgs-bd]

@jlarsen-usgs many thanks!

[dbrakenhoff]

With the update in PR #2285, the resample filters out NaN values prior to resampling using nearest interpolation.

A user asked me whether this is always the intended behavior? They wanted to resample a raster with layer elevations (including NaNs where this particular geological feature is not present). The result should be layer elevations on a finer grid with NaNs where the original dataset contained NaNs. The NaNs in the result can then be used to set the idomain for that layer. The previous version of flopy returned a result containing NaNs, though that might not have been the intended behavior judging from the comments above.

So my question:

• How would one resample to a finer grid, and keep the NaNs from the original dataset? Or is this not possible with the current methods in flopy?

[jlarsen-usgs]

@dbrakenhoff

The update to the "nan" filtering behavior applies to the interpolation based resampling (bi-linear, bi-cubic, and nearest neighbor).

The other geostatistical methods, such as min, max, median, mean, and mode do not filter out the raster nodata values. My general recommendation for working with DEM data is to sample based the minimum elevation. When the grid discretiztion is smaller than the raster discretization it will set grid cells to the nearest neighbor with any these methods, because grid cell/raster cell intersection would produce a single value. Here is the test code I used to confirm this

from flopy.discretization import StructuredGrid
import numpy as np
from flopy.utils import Raster
from flopy.plot import PlotMapView
from pathlib import Path
import matplotlib.pyplot as plt
from shapely.errors import ShapelyDeprecationWarning
import warnings
warnings.simplefilter("ignore", ShapelyDeprecationWarning)


def create_grid(raster, dx, dy, buffer=0):
    nlay = 1
    xmin, xmax, ymin, ymax = raster.bounds
    adj = buffer / 2
    xmin -= buffer
    ymin -= buffer
    xmax += buffer
    ymax += buffer

    nrow = int(np.ceil((ymax - ymin) / dy))
    ncol = int(np.ceil((xmax - xmin) / dx))

    delc = np.full((nrow,), dy)
    delr = np.full((ncol,), dx)
    top = np.ones((nrow, ncol))
    botm = np.zeros((nlay, nrow, ncol))
    idomain = np.ones(botm.shape, dtype=int)

    return StructuredGrid(
        delc=delc,
        delr=delr,
        top=top,
        botm=botm,
        idomain=idomain,
        xoff=xmin,
        yoff=ymin,
    )


if __name__ == "__main__":
    # from the flopy repository base directory
    dem_file = Path("./examples/data/options/dem/dem.img")
    rstr = Raster.load(dem_file)

    test_grid = create_grid(rstr, 15, 15, buffer=1000)
    resampled_dem = rstr.resample_to_grid(
        test_grid,
        band=rstr.bands[0],
        method="min",
        extrapolate_edges=False
    )

    fig, (ax0, ax1) = plt.subplots(ncols=2)

    rstr.plot(ax=ax0)
    pmv0 = PlotMapView(modelgrid=test_grid, ax=ax0)
    pmv0.plot_grid(lw=0.1, alpha=0.7)

    pmv1 = PlotMapView(modelgrid=test_grid, ax=ax1)
    pmv1.plot_array(resampled_dem, masked_values=rstr.nodatavals)
    pmv1.plot_grid(lw=0.1, alpha=0.7)

    plt.show()

[dbrakenhoff]

Thanks @jlarsen-usgs, I'll check this out!

[Arjanvw]

Hi @jlarsen-usgs ,

I'm the user dbrakenhoff is helping. The solution provided seems to work with upscaling. We want to use Nearest Neighbor since we want to resample from 100x100m to 20x20m, so downscaling. See the example script below (based on your example). Method "min" only provides results on the cell centers, method "nearest" provides a filled raster. However, I need the NaN values to set the idomain of the clay layer to inactive. Is it possible to add an option "NaN-in=NaN-out" for the "nearest", "linear" and "cubic" methods?

# %% Load python packages

import numpy as np
import matplotlib.pyplot as plt

from flopy.plot import PlotMapView
from flopy.utils import Raster
from flopy.discretization import StructuredGrid

import flopy

print(f"flopy version: {flopy.__version__}")

# %% ascii file for clay layer

with open('top_clay.asc', 'w') as f:
    f.write('ncols 4\n')
    f.write('nrows 4\n')
    f.write('xllcorner 0.0\n')
    f.write('yllcorner 0.0\n')
    f.write('cellsize 100\n')
    f.write('nodata_value -9999\n')
    f.write('0.1 0.2 0.3 0.4\n')
    f.write('0.2 0.3 -9999 0.5\n')
    f.write('0.3 0.4 -9999 0.6\n')
    f.write('0.4 0.5 0.6 0.7\n')

rstr = Raster.load("top_clay.asc")

# %% make dis

dx = 20
dy = 20

nlay = 1
xmin, xmax, ymin, ymax = rstr.bounds

nrow = int(np.ceil((ymax - ymin) / dy))
ncol = int(np.ceil((xmax - xmin) / dx))

delc = np.full((nrow,), dy)
delr = np.full((ncol,), dx)
top = np.ones((nrow, ncol))
botm = np.zeros((nlay, nrow, ncol))
idomain = np.ones(botm.shape, dtype=int)

test_grid = StructuredGrid(
    delc=delc,
    delr=delr,
    top=top,
    botm=botm,
    idomain=idomain,
    xoff=xmin,
    yoff=ymin,
)

# %% test resample methods

resampled_min = rstr.resample_to_grid(
    test_grid,
    band=rstr.bands[0],
    method="min",
    extrapolate_edges=False
)

resampled_nn = rstr.resample_to_grid(
    test_grid,
    band=rstr.bands[0],
    method="nearest",
    extrapolate_edges=True
)

# %% plot
fig, (ax0, ax1, ax2) = plt.subplots(ncols=3,figsize=(15, 5))

rstr.plot(ax=ax0)
pmv0 = PlotMapView(modelgrid=test_grid, ax=ax0)
pmv0.plot_grid(lw=0.1, alpha=0.7)
ax0.set_title('input raster')

pmv1 = PlotMapView(modelgrid=test_grid, ax=ax1)
pmv1.plot_array(resampled_min, masked_values=rstr.nodatavals)
pmv1.plot_grid(lw=0.1, alpha=0.7)
ax1.set_title('method min')

pmv2 = PlotMapView(modelgrid=test_grid, ax=ax2)
pmv2.plot_array(resampled_nn, masked_values=rstr.nodatavals)
pmv2.plot_grid(lw=0.1, alpha=0.7)
ax2.set_title('method nearest')

plt.tight_layout()
plt.show()

[vincentpost]

The problem with 'min' only providing results at the cell centers of the raster data can be avoided by calling rasterstats.zonal_stats with all_touched=True. By default it is set to False, see https://pythonhosted.org/rasterstats/manual.html#rasterization-strategy for the different behaviors. I am not sure if there is any harm in setting it to True by default? Or perhaps it can become a kwarg for resample_to_grid.

I agree with @Arjanvw that it would be good to also have the choice to either keep or remove the nodatavals from the data before interpolation with griddata. Looking at the code to see how it could be implemented, I was wondering though why nodatavals are set to nan for 'nearest' and 'linear' but not for 'cubic'. I was afraid to break it so I didn't touch it.

[jlarsen-usgs]

@Arjanvw, @dbrakenhoff

I just opened a PR (#2331) to revert the behavior of resample_to_grid() back, as there were some unintended consequences that occured with the previous PR that addressed filling nan values. With this PR, the extrapolate_edges kwarg now works with the nearest neighbor method and can be optionally used to fill nan values. If extrapolate_edges=False (default) nan values are maintained, which was the behavior prior to #2285.