pointrasterterrset
I know it is possible to use tools like ArcGis's 'extract sample' tool to extract values at a given number of points. I would like to sample every pixel in a raster, and drill down through all other rasters (which have the same extent and resolution). How could I do this? I can use QGIS, Arc, or Idrisi.
Would the best way be to generate regular points at a spacing of 'RESOLUTION' and an initial offset of 'RESOLUTION/2'?
Using R, once you have read in or created rasters x and y, the regression is performed and a raster of residuals r can be produced with three commands:
model <- lm(getValues(y) ~ getValues(x), na.action=na.exclude)
r <- y
r[] <- residuals.lm(model)
(To verify everything is working correctly, printing and plotting the results is a good idea.) The workflow generalizes readily to multiple independent rasters and to other models such as generalized linear models, nonlinear fits, models with spatial correlation, etc.
Let's create some sample data:
#
# Create raster objects.
#
nRows <- 7
nCols <- 10
x <- raster(nrow=nRows, ncol=nCols)
y <- x
e <- x
#
# Simulate a bivariate relation y ~ x with error e.
#
set.seed(17)
e[] <- rnorm(nRows*nCols)
x[] <- 1:nCols
y <- x * 0.5 + 5.0 + e
Thus, the independent raster is x and the dependent is y; the errors in the linear relation are e. (In practice the errors are never known--they have to be estimated by the residuals--but with a simulated dataset we can hold on to the error values for checking things later.)
Do the regression:
model <- lm(getValues(y) ~ getValues(x), na.action=na.exclude)
summary(model)
The output lets us check that everything is working as expected. Here's part of it:
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.92223 0.25842 19.05 <2e-16 ***
getValues(x) 0.53734 0.04165 12.90 <2e-16 ***
...
Residual standard error: 1.001 on 68 degrees of freedom
The estimates are highly significant and quite close to the actual values (4.92 vs. 5.00 and 0.54 vs. 0.05). The residual standard error of 1.001 closely estimates the error built into the simulation, 1.0. Finally, 68 = 7*10 - 2 counts the cells minus the two estimated parameters. Everything looks good.
Let's create your grid of residuals and plot it:
r <- y # Create a grid object for the residuals
r[] <- residuals.lm(model)
plot(r)
As a double-check (because we have been copying values out of and then back into raster format), let's plot the residuals against the true errors. The points should closely follow the line y=x, with just a tiny bit of scatter:
plot(getValues(e), getValues(r), xlab="Error grid", ylab="Residual grid")
Everything is fine.
The simple way to do this in QGIS is to use the Raster Calculator (Raster->Raster Calculator). You have a couple of options. The easiest to explain/understand is to make a unitary raster from your mask (all data set to either 1 or NoData) and then multiply your clip layer by the unitary mask layer.
To ensure the extents match the mask layer, in the raster calculator window select the mask layer in the 'Raster bands' list on the left and then click the 'Current layer extent' button on the right.
You can create a unitary mask on the fly by using a conditional statement (see the link) something like this:
(maskLayer@1 >= 0) * clipLayer@1
This statement basically says: treat everything in my mask layer that is not NoData as being equal to 1 (NoData stays as NoData). Just be sure to remember to set the extent (see above)!
Raster Calculator and paletted data:
Any output from the raster calculator will be just values and not carry over any information contained within a colour pallet. You have a couple of options to 'get your palette back':
Best Answer
The Raster to Point tool in ArcGIS will create a point at the centre of each cell. Then you can use Sample.
You could change the resolution of your raster and create points for that with RESOLUTION/2 and combine the points. But I don't see the need. Any point created with resolution/2 will fall within the original pixel. Or on the edge of two pixels.