Focal regression using python

focal_regression.py

"""
This code uses the rasterio library to open the two input rasters and the scipy library
to perform the focal regression using the generic_filter function. The size parameter 
specifies the size of the neighborhood around each pixel to use for the regression, 
in this case a 3x3 window. The extra_arguments parameter is used to pass the second raster 
to the np.polyfit function, which fits a polynomial of degree 1 (a linear regression) to 
the two rasters. The result is a 3D array with the slope and intercept values for each pixel. 
The code then extracts the slope and intercept arrays from the result and saves them to file 
using rasterio.

Note that this code assumes that the two rasters have the same spatial reference and dimensions. 
You may need to modify the code to handle rasters with different spatial references or dimensions.
"""
import rasterio
import numpy as np
from scipy import ndimage

# Open the two rasters
with rasterio.open('raster1.tif') as src1:
    raster1 = src1.read(1)

with rasterio.open('raster2.tif') as src2:
    raster2 = src2.read(1)

# Perform focal regression
result = ndimage.generic_filter(raster1, np.polyfit, size=(3,3), extra_arguments=(raster2, 1))

# Extract the slope and intercept from the result
slope = result[:,:,0]
intercept = result[:,:,1]

# Save the slope and intercept rasters to file
with rasterio.open('slope.tif', 'w', **src1.meta) as dst:
    dst.write(slope, 1)

with rasterio.open('intercept.tif', 'w', **src1.meta) as dst:
    dst.write(intercept, 1)