I have two rasters of the same region and of same resolution and I want to plot in a scatter plot their values and compute the best - fit line with r-square. Any idea how can I do this on Arcmap or QGIS or extract the values of the rasters in a file so I can import it into Excel? Unfortunately I don't use R neither python.
You can use the Image Correlation, Image Regression, and Feature Space Plot tools in Whitebox Geospatial Analysis Tools to achieve this. It works on whole images, with the caveat that when you have very high sample sizes (i.e. millions of pixels) even very small differences will yield statistical significance, which may not be physically meaningful.
Here's an example of an Image Regression output for two bands of Landsat 8 imagery:
Doing regression entirely within a GIS is a bit of a tour de force, but it is possible (and perhaps has some merit when incorporated within a longer workflow that would otherwise be interrupted by the use of a statistical computing platform.)
Whenever you have a huge set of matched data (X,Y), such as afforded by two collocated rasters, then the ordinary least-squares fit is best found with a series of simple steps (which avoid floating point errors that can occur when applying conventional formulas):
Find the averages of the X and Y. Call these mX and mY, respectively. (Ordinarily these would be represented as rasters having constant values in their cells.)
Recenter X and Y by subtracting their averages from each. (That is, change X to X-mX and Y to Y-mY).
Find the averages of the squared recentered values: these are the variances of the rasters. (Again, these would be represented as constant rasters.)
Divide the recentered values by the square roots of their variances. (Call these square roots sX and sY, respectively.) The resulting rasters have been standardized. Most of their values will lie between -2 and 2, with a few extending beyond these limits (sometimes far beyond).
Multiply the two standardized rasters. The average of these products is the "beta" for the regression: it is the standard coefficient. Its square is R-squared. Because this is the desired answer, if you won't be using it for prediction or other processing, you might prefer to obtain it in a (one-line) table rather than representing it as a raster.
Incidentally, the slope of the regression line of Y against X is obtained by converting back from the standardized to the original values: this will multiply Y by sY and X by sX. Thus, beta has to be multiplied by sY and divided by sX. The slope of the regression line of X against Y similarly is obtained by dividing beta by sY and multiplying it by sX. (Those two slopes will not equal each other unless beta is at the extreme value of 1 or -1.)
Only two kinds of raster operations are involved in all this: zonal statistics (to average values across a raster) and local arithmetic (to subtract values from each cell, square the values, take square roots, multiply values in two grids, and so on). To carry out the zonal statistics, first create a zone grid by overlaying the X and Y rasters and setting the resulting values to a constant (such as 0 or 1); every cell without a value in both X and Y will be null. This makes the entire raster into a single zone over which all averages will take place.
To draw the scatterplot, consider sampling from the rasters in order to limit the points to a manageable number.
Usually you'll have enough information for such a calculation using a few thousand of points, and excel might be limited if you use the full dataset.
So you could generate a set of random points (simple random sampling or regular grid with a random origin), extract the pixel values for each point and then export the attribute table to excel for further calculation.