Updating all elements with a given value in a 2darray with neighboring values average

Viewed 114

Consider a 2darray numeric array where a specific value, i.e. 0, represents "no data". In some instances, it may be desirable to update these no data values with the neighboring cell's value. For example, I have two grids laid on top of each other - one has parcel IDs and the other has transportation times to major cities. Some IDs will only overlap with no data transportation time cells. Here, it would actually be an ok guess to just use the average of neighboring grid cells. Focusing more on the abstract question, though, consider the following set up

v_shp_ids = np.array([[1, 1, 1, 1],
                      [1, 1, 1, 1],
                      [1, 1, 1, 1],
                      [1, 1, 1, 1]])

v_shp_ids = v_shp_ids.ravel()

raster = np.array([[0, 0, 0, 0],
                   [0, 0, 0, 0],
                   [0, 0, 1, 0],
                   [0, 0, 0, 0]])

raster_r = raster.ravel()

zs = pd.concat([pd.Series(v_shp_ids), pd.Series(raster_r)], axis=1)

Here, v_shp_ids is saying Parcel ID 1 is represented by all the cells. raster is saying there is only one cell with a value. I want to get raster filled with 1's. (Starting with a simple case).

This is my solution

 # Dictionary to store no data updates
na_update = {}

# Num rows in raster 2darray
i_end = len(raster)
# Num cols in raster 2darray
j_end = len(raster[0])
# Hard coding nodata = 0 for the example arrays I input above
nodata = 0

# Number of raster cells with no data
fill_length = len(zs.loc[zs[1] == nodata])
# Copy the dataframe for no data raster cells only
zs_na = zs.loc[zs[1] == nodata].copy()

# While the updated data dictionary is not full, loop
while len(na_update) != fill_length:
    # For each 1darray index in the df
    for ind, row in zs_na.iterrows():
        # get the corresponding raster grid cell
        i = ind // j_end
        j = ind % j_end
        # get the slice of neighboring pixel cells and flatten
        neighb = raster[max(i-1,0):min(i+2,i_end),
                        max(j-1,0):min(j+2,j_end)].flatten()
        # get the average of these neighboring, excluding na values
        neighb_avg = np.ma.masked_values(neighb, nodata).mean()
        # store average in dictionary if we can and update the pixels
        if neighb_avg is not np.ma.masked:
            na_update[ind] = neighb_avg.astype(raster_r.dtype)
            raster[i, j] = na_update[ind]
    # update the df to the remaining pixels with na remaining
    zs_na = zs_na[~zs_na.index.isin(na_update.keys())]

I'm aware that convolution could be an approach here but some of my concerns are a) I don't want to calculate nearest neighbor averages for every cell since in many instances I'm not dealing with no data cells b) I can't ex-ante tell what neighborhood size to use to ensure all no data cells get a value and to ensure they all do, the neighborhood might be too big to continue being a reasonable approximation (Sure, you might think the whole method is not a reasonable approximation)

I have two main questions

  1. Is there a way to use indexing on the dataframe (or even the raster grid) to do the same thing? I had played around with variations of zs.loc[zs[1] == nodata, 1] = # Trying to take the average of adjacent elements in the raster grid

  2. Have I completely missed a reasonable method for filling no data cells?

I found this a fun problem to work on and I'm very much looking forward to ideas/discussion.

1 Answers

Here is a rasterio function that works great: https://rasterio.readthedocs.io/en/latest/api/rasterio.fill.html

The function is fillnodata from the rasterio.fill module. From the description:

Fill holes in raster data by interpolation

This algorithm will interpolate values for all designated nodata pixels (marked by zeros in mask). For each pixel a four direction conic search is done to find values to interpolate from (using inverse distance weighting). Once all values are interpolated, zero or more smoothing iterations (3x3 average filters on interpolated pixels) are applied to smooth out artifacts.

The function call is

rasterio.fill.fillnodata(image, mask=None, max_search_distance=100.0, smoothing_iterations=0)

And I passed in

rasterio.fill.fillnodata(image = np.ma.masked_array(raster, nodata))
Related