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.
3 Answers
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:
-
+1. Note, however, that the p-value will be meaningless with many rasters (such as those created via interpolation or density calculations) because they cannot adequately be modeled as collections of variables with independent errors. Much more effective for rasters would be a dynamic windowing tool that displays a measure of relationship within local (smaller) neighborhoods, allowing the viewer to explore how their relationship varies with location.– whuberCommented Aug 15, 2014 at 13:47
-
Hi @whuber, yes, that is true. I'm very interested in your idea of the dynamic windowing tool for local p values. Do you have a reference for this? Thanks. Commented Aug 15, 2014 at 14:32
-
I have no reference. This has echoes of GWR, but it's a little different. You could carry it out by following the recipe in my reply, replacing all zonal statistics with focal statistics (which would be quickly computed using the FFT since they are all convolutions). The result would be a grid of beta values, or of R-squareds, or of slopes, or anything else you care to derive from them. I wouldn't bother to compute p-values--they would be practically meaningless--but the resulting raster could be an extremely efficient and informative way to explore nonstationary relationships.– whuberCommented Aug 19, 2014 at 22:15
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.
-
this is an excellent approach. My only addition is that instead of sampling the raster to make the scattergram, create a 2D histogram, e.g. oceanpython.org/2013/02/25/2d-histogram Commented Aug 14, 2014 at 23:52
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.
In ArcGIS, you'll need either "fishnet" or "generate random points" for the first step, then "extract multiple value to point" (spatial analyst licence required) and finally "table to table"
-
@radouxjo thank you for your solution, but I really need to compare the whole raster with the other raster... Commented Aug 14, 2014 at 21:01
-
statistically this will not change much, and this is typically something that I would not do with ArcGIS and Excel, but if you really need it you can also convert one of the raster to point (using the tool of the same name) as the first step, then follow with the second and third step.– radouxjuCommented Aug 14, 2014 at 21:06