Linear regression between every 3×3 pixels of two rasters

I have two rasters of land surface temperature and elevation, LST and ELEV. They have equal number of rows and columns (same boundary, same number of pixels).I want to create a new raster of regression slope or temperature lapse rate (RSLOPE) with same number of rows and columns as LST and ELEV. I want to create this raster by performing linear regression between every 3×3 pixels of the two rasters (LST and ELEV) such that each pixel of the RSLOPE will hold the regression slope value obtained from linear regression of the corresponding 3×3 pixels in LST and elev that surround that pixel.

Answer 1

[Guillaume]

For what you want to do, blockproc wouldn't work as it does not overlap the blocks. nlfilter and colfilt do on the hand, but none of these functions take two images as inputs. In addition nlfilter and colfilt both do zero-padding at the end which is not what you want.

What I don't understand is why you only expect one image out of your linear regression. You usually get two terms out of a linear regression, the slope and the offset. The equation of the line in your image is incorrect, it's actually LST = RSLOPE * ELEV + ROFFSET.

Regardless, what you could use is im2col to convert the blocks of both input images into columns, perform your regression on each column, then reshape the result back to an image.

collst = im2col(LST, [3 3]);
colelev = im2col(ELEV, [3 3]);
rslope = zeros(1, size(collst, 2));
roffset = rslope;
%perform linear regression on each column
for col = 1:size(collst, 2)
  p = polyfit(colelev(:, col), collst(:, col), 1);
  rslope(col) = p(1);
  roffset(col) = p(2);
end
%reshape result  
rslope = reshape(rslope, size(LST, 1) - 2, size(LST, 2) - 2);
roffset = reshape(roffset, size(LST, 1) - 2, size(LST, 2) - 2);

The slowest part of the above is going to be the polyfit applied to each column one by one. As it's a linear regression, the equations for the slope and offset are well know so you could implement them directly instead of using polyfit. The advantage of that would be much faster processing as you could apply the equation directly to the whole collst and colelev matrices without a loop.

You'll note that the result is 2 rows and columns shorter than the inputs. While im2col doesn't do padding, it also does not do shrinking of the blocs, so the edges are ignored. I don't have a solution for that.

[Alejandro]

I am sorry for my incomplete figure. I do understand that I would get two terms out of a linear regression, the slope and the offset. I was trying to highlight that I am only interested in the RSLOPE value and then converting those RSLOPE values obtained after linear regression of pixel pair in LST-ELEV pair into a raster of RSLOPE values with pixel values corresponding to the pixels in LST-ELEV pair which were involved in linear regression. I also need the RSLOPE raster exactly equal to the size of LST or ELEV raster. Thank you so much for your codes and your suggestions again. I wonder if 'reshape result' would exactly place an RSLOPE pixel in the same order (same row,same column) of the LST-ELEV pixel pair involved in the linear regression to obtain that particular RSLOPE pixel. Further, I wonder if I may take the two rasters PPT and ELEV occupying a larger area than my area of interest and apply what you have suggested to calculate the regression slope and obtain the raster of regression slope (RSLOPE) and later clip this raster RSLoPE using shapefile of my area of interest to avoid the limitations of 'im2col' (limitations: result is 2 rows and columns shorter, edges are ignored). Thanks again for your prompt reply and your kind consideration. 

[Guillaume]

I wonder if 'reshape result' would exactly place an RSLOPE pixel in the same order (same row,same column) of the LST-ELEV pixel pair involved in the linear regression

The answer to that is yes, and it is guaranteed by the documentation. See the tips section of im2col.

I wonder if I may take the two rasters PPT and ELEV occupying a larger area than my area of interest

That would be the easiest way to obtain a result of the right size. The easiest way to do this is to padarray the inputs with a padsize of [1 1]. However, whichever constant you'd pad with would affect the regression.

One option would be to pad with NaN and remove these NaN before the polyfit call:

collst = im2col(padarray(LST, [1 1], NaN), [3 3]);
colelev = im2col(padarray(ELEV, [1 1], NaN), [3 3]);
rslope = zeros(1, size(collst, 2));
%perform linear regression on each column
for col = 1:size(collst, 2)
  x = colelev(:, col);
  y = colelev(:, col);
  keep = ~isnan(x);
  p = polyfit(x(keep), y(keep), 1);
  rslope(col) = p(1);
end
rslope = reshape(roffset, size(LST));

However, it will be much more difficult to vectorise the regression if there are NaNs present.

[Guillaume]

Actually, thinking more about it, adding the NaN wouldn't make the vectorisation of the regression more complicated since sum, mean etc. now all have an 'omitnan' option.

Completely untested, using the formula from here:

collst = im2col(padarray(LST, [1 1], NaN), [3 3]);
colelev = im2col(padarray(ELEV, [1 1], NaN), [3 3]);

%perform linear regression
sumxy = sum(collst .* colelev, 'omitnan');
sumx = sum(colelev, 'omitnan');
sumy = sum(collst, 'omitnan');
sumxsq = sum(colelev.^2, 'omitnan');
n = sum(~isnan(collst));
slope = (n .* sumxy - sumx .* sumy) ./ (n .* sumxsq - sumx.^2);  

%reshape
slope = reshape(slope, size(LST));

This is a lot faster than polyfit over each column.

[Alejandro]

Forgive me for taking such a long time to reply. I took a little time to figure out the codes written by you. Many thanks for your formula. It works absolutely fine. This was really very helpful. This thread will be very useful to someone trying to work with low resolution climate data in high mountains where data is scarce. I still have a query. Can we write the 'slope' matrix created at the end into a raster (probably tiff file)? I want to perform spatial analysis with this data. Thanks again for your kind support.

[Guillaume]

"This thread will be very useful to someone ..."

Then I suggest you accept the answer.

"Can we write the 'slope' matrix created at the end into a raster"

The easiest is to rescale it to (1, 255) and save it as indexed

rasterslope = 1 + (slope - min(slope(:)) / (max(slope(:)) - min(slope(:)) * 255;
imwrite(rasterslope, parula(256), 'somename.png');

[Alejandro]

I am sorry I did not realise that I have not already accepted the answer. I converted the matrix into ASCII file and then converted ASCII to raster using ArcGIS and then later georeferenced it. Many thanks again for your kind support and serious consideration.This was really helpful and needed.