How to select one band in image using rasterio

Question

I have two different images and I want to select one specifi band from each one of them and calculate regression between them. Each one of those bands is a result of calculation I did (in the first image , band' no. 4 is the results of NDVI calculation, and the other band is SAR db values). The problem is that I haven't found a way to do it.

I used "read" in order to try to select one band but it didn't work when I tried to do things with this.

My end goal is to actually calculate linear regression between those two bands ,which I also still don't know how to do in rasterio library.

my code:

import pandas as pd
import rasterio
import numpy as np
from rasterio.plot import show
import seaborn as sns
%matplotlib inline

#NDVI image

NDVI = rasterio.open(r"C:\Users\PATH\TO\IMAGE\python\NDVI.tif")

show(NDVI)

#SAR image 

SAR = rasterio.open(r"C:\Users\path\to\IMAGE\python\SAR.tif")

NDVIband=NDVI.read(4)
SARband=SAR.read(2)

#from here nothing work and I feel lost. 
#I thought to use this library but seems like it works only with pandas (?)
from sklearn.linear_model import LinearRegression

as i'm very new to python and image processing, any clue will be valuable for me ...

Answer 1

I don't have a SAR image so, I assumed that this image has same values as NDVI image. Result of this kind of linear regression, approximately, would be a straight line with slope 1 and intercept 0. You don't have any result because you first need to prepare numpy arrays by using flatten and reshape numpy methods. My code, by using this link as guide for scikit-learn python module is:

import pandas as pd
import rasterio
import numpy as np
from rasterio.plot import show

#NDVI image (independent variable) 
NDVI = rasterio.open("/home/zeito/pyqgis_data/NDVI.tif")
#show(NDVI)

#SAR image (assuming NDVI = SAR) as dependent variable 
SAR = rasterio.open("/home/zeito/pyqgis_data/NDVI.tif")

NDVIband = NDVI.read(1).flatten().reshape((-1, 1))
SARband = SAR.read(1).flatten()

from sklearn.linear_model import LinearRegression

model = LinearRegression().fit(NDVIband, SARband)

r_sq = model.score(NDVIband, SARband)

print('coefficient of determination:', r_sq)
print('intercept:', model.intercept_)
print('slope:', model.coef_)

After running above code in Python Console, as expected, result was:

coefficient of determination: 0.9999999999992975
intercept: 2.682209e-07
slope: [0.99999917]