Outlier's neighbours of a raster analysis in Python/GIS

Question

Starting with point clouds from various drone flights, I created a 0'1x0'1m cell raster with the DEM of Differences (DoD). Next, in order to clean outliers cells, I calculated the 0'1 and 99'9 percentile with Python:

from osgeo import gdal, gdal_array
import numpy as np

raster = r"E:\[...]\DoD.tif"

nanvalue = -3.4028230607370965e+38 #TIF's NaN value

perc = [0.1,99.9] #Percentiles

rArray = gdal_array.LoadFile(raster) #Read raster as numpy array

mask = np.ma.masked_where(rArray == nanvalue, rArray) #Mask for NaN
mask = np.ma.filled(mask, np.nan)

for p in perc:
    print('{0}th percentile: {1}'.format(p, np.nanpercentile(mask,p)))

Knowing these values, I need to study the array neighbours of each out-of-percentile cell, so I see if really is an outlier or could be real peak value. If the cell is isolated, that cell will be deleted but, if the cell has 3 or more neighbours with a similar value (i.e. standard deviation of the DoD), the cell remains.

I could list every value out of range with a loop/conditional, but I'm struggling with making it analyze its neighbours. Some ideas?

Also, I'm doing this in Python to make the process repeatable, but I'm open to use QGIS/ArcGIS tools too if you know any.

Answer 1

Well, I've been trying to solve the problem. For now, the best aproximation I could have was this one.

Basically, mask the NaNs in the first array, calculate the percentiles and the bounds (percentile value +- sd). All outliers will be set as TRUE. Then, amplify the TRUE values around the outlier cells (3x3 matrix), so I can compare the neighbors’ value to the outlier value. If the difference is in the +-sd bound in 2 or less neighbors, the outlier will be considered as isolated pixel and will be deleted.

I’m going to post it because I have seen few questions and answers about it. Hope it helps someone.

from osgeo import gdal, gdal_array
import numpy as np
from scipy.ndimage import label, binary_dilation
import matplotlib.pyplot as plt

raster = r"E:\[...]\DoD.tif"

#TIF's NaN value
nanvalue = -3.4028230607370965e+38 

#Percentiles
perc = [0.1, 99.9] 

#Read raster as numpy array
rArray = gdal_array.LoadFile(raster) 

#Mask for NaN
mask = np.ma.masked_where(rArray == nanvalue, rArray) 
mask = np.ma.filled(mask, np.nan)

for p in perc:
    print('{0}th percentile: {1}'.format(p, np.nanpercentile(mask, p)))

#####Outliers and neighbors#####

#Outliers lower and higher than percentiles along SD threshold
lower_bound = np.nanpercentile(mask, perc[0])
upper_bound = np.nanpercentile(mask, perc[1])
sd_threshold = np.nanstd(mask)

outliers = np.logical_or(mask < lower_bound - sd_threshold, mask > upper_bound + sd_threshold)

#Print the bounds
print("Standard Deviation Threshold:", sd_threshold)
print("Lower bound - SD:", lower_bound - sd_threshold)
print("Upper bound + SD:", upper_bound + sd_threshold)

#Turn all outlier's neigbors to TRUE in a 3x3 matrix
neighbors_array = binary_dilation(outliers, structure=np.ones((3, 3)))

#TRUE checks. Not all outliers are completely surrounded, so could be less than 8 times the outliers
num_outliers = np.sum(outliers)
num_neighbors = np.sum(neighbors_array)
print("Number of TRUE elements in outliers:", num_outliers)
print("Number of TRUE elements in neighbors_array:", num_neighbors)

#Remove outliers with 2 or less neighbors with similar value (+-sd_threshold)
labeled_array, num_labels = label(neighbors_array) #To label elements

def clean_outliers_by_nb(label_id):
    #Label every connected element. If label_id = 1, all component's ones will be TRUE
    component = (labeled_array == label_id) 
    neighbors_values = mask[component]
    
    #Count neighbors with the outlier value within the threshold
    num_valid_neighbors = np.sum(np.abs(neighbors_values - mask[component]) <= sd_threshold)
    
    #If 2 or lower valid neighbors, remove the outlier
    if num_valid_neighbors <= 2:
        outliers[component] = False
        
#Loop the function
for label_id in range(1, num_labels + 1):
    clean_outliers_by_nb(label_id)

#TRUE check
num_outliers = np.sum(outliers)
print("Number of TRUE elements in outliers after removal:", num_outliers)

#####Save as raster and txt#####

#Copy of the original array to store and apply changes to the modified array
output_raster_path = r"E:\[...]\DoD_clean.tif"
modified_array = mask.copy() 
modified_array[outliers] = np.nan #Outliers as NaN 

#Original raster's info
geotransform = gdal.Open(raster).GetGeoTransform()
projection = gdal.Open(raster).GetProjection()

#Create a new raster with the modified array
driver = gdal.GetDriverByName("GTiff")
modified_raster = driver.Create(output_raster_path, rArray.shape[1], rArray.shape[0], 1, gdal.GDT_Float32)

#Set info
modified_raster.SetGeoTransform(geotransform)
modified_raster.SetProjection(projection)

#Write the modified array to new raster
modified_raster.GetRasterBand(1).WriteArray(modified_array)

modified_raster = None #Close raster dataset

#Save txt
output_txt_path = r"E:\[...]\DoD_clean.txt"
np.savetxt(output_txt_path, modified_array, fmt='%.8f', delimiter=' ')

#####Plot#####

modified_raster = gdal.Open(output_raster_path)
modified_array = modified_raster.GetRasterBand(1).ReadAsArray()

fig = plt.figure(figsize=(8, 6), dpi=300)
plt.imshow(modified_array, cmap='turbo_r')
plt.axis('off')
plt.colorbar(label='m')
plt.title('DEM of Difference')
plt.show()