This lab will help you to estimate the uncertainty in your data.
1.Obtain a continuous raster. This can be one from your own research or the Humboldt County 10 meter DEM.
2. Use the code below to create a function that will compute the mean or standard deviation of a down-sampled raster. Note that the code contains examples of if statements, while loops and functions in R.
The code below requires the packages rgdal and raster. Please see the R website for information on installing these packages. I also highly recommend using small rasters to get started as these functions will take much longer with large rasters.
###############################################################################
# Function to obtain a down-sampled raster where each pixel represents the
# mean of the pixels that it overlapped with in the original raster.
#
# Inputs:
# - OriginalRaster - The raster to be downsampled
# CellSize - The number of pixels in the x and y directions that will be
# downsampled for each new pixel(i.e. a 2 to convert 4 pixels to 1)
# Type - The type of raster to be created
# MEAN for a raster with the average value for the sub-sampled pixels
# STDDEV for a raster with the standard deviation for the sub-sampled pixels
# MIN for a raster with the minimum values from the sub-sampled pixels
# MAX for a raster with max values
# FIRST for the first value in each sub-sample.
###############################################################################
###############################################################################
# Constant definitions
###############################################################################
# these constants specify the "Type" of downsampled raster to obtain
MEAN=1
STDDEV=2
MIN=3
MAX=4
FIRST=5
###############################################################################
# The Function
###############################################################################
DownSampleRaster=function(OriginalRaster,CellSize,Type)
{
# Get the matrix containing the pixels for the original raster
TheMatrix=as.matrix(OriginalRaster)
# Get the number of rows and columns in the original raster
NumRasterColumns=ncol(OriginalRaster)
NumRasterRows=nrow(OriginalRaster)
# Get the number of rows and columns in the new raster
NumRows=NumRasterRows/CellSize
NumColumns=NumRasterColumns/CellSize
NumRows=as.integer(NumRows)
NumColumns=as.integer(NumColumns)
# Create the matrix that will include the pixels for the new raster
NewMatrix=matrix(nrow = NumRows,ncol=NumColumns) #
# Setup a variable with the number of cells in the subset as we'll use this again and again
NumSubsetPixels=CellSize^2
# Initialize the row and column counters
Row=0
Column=0
while (Row<NumRows) # interate through each row in the new raster (we skip the last row if partial)
{
Column=0 # reset the column counter to go through each column for each row
while (Column<NumColumns)
{
# Find the row and column indexes into the original raster
FirstPixelX=Column*CellSize+1 # R indexes from 1
FirstPixelY=Row*CellSize+1
LastPixelX=FirstPixelX+CellSize-1
LastPixelY=FirstPixelY+CellSize-1
# Get the Subset from the original raster that is downsampled into the new raster
Subset=TheMatrix[FirstPixelY:LastPixelY,FirstPixelX:LastPixelX]
# Find the mean value
# This can be replaced by other code to find the minimum, maximum, and std. dev.
if (Type==MEAN)
{
Value=mean(Subset)
}
else if (Type==STDDEV)
{
Value=sd(Subset)
}
else if (Type==MIN)
{
Value=min(Subset)
}
else if (Type==MAX)
{
Value=max(Subset)
}
else if (Type==FIRST)
{
Value=sd(Subset[1])
}
# Save the new value into the raster
NewMatrix[Row+1,Column+1]=Value
# Move to the next column
Column=Column+1
}
Row=Row+1
}
# Create the new raster
NewRaster=raster(NewMatrix)
# Copy the spatial reference from the original raster to the new raster
TheProjection=projection(OriginalRaster)
projection(NewRaster)=TheProjection
# Copy the raster extent (even downsampled, the extent is the same)
# The pixels do get larger
xmin(NewRaster)=xmin(OriginalRaster)
xmax(NewRaster)=xmax(OriginalRaster)
ymin(NewRaster)=ymin(OriginalRaster)
ymax(NewRaster)=ymax(OriginalRaster)
# return the new raster
return(NewRaster)
}
3. Use the code below to call the function above to complete a mean or standard deviation raster by passing a defined "TYPE" to in with the function. In a standard deviation raster, each pixel represents the standard deviation of the pixels down sampled from in the original pixel. Try this for at least 3 cell sizes.
# Setup the folder and file paths FolderPath="C:\\jg2345\\GSP 570\\Lab 3 Predictor Uncertainty\\" #InputFilePath=paste(FolderPath,"HumboldtDEM_10Meter.tif") InputFilePath=paste(FolderPath,"Cropped_1.tif",sep="") # Read the raster into a raster object OriginalRaster=raster(InputFilePath) # Setup the cell size CellSize=2 # Get the downsampled raster NewRaster=DownSampleRaster(OriginalRaster,CellSize,STDDEV) # Setup the output file path OutputFilePath=paste(FolderPath,"NewRaster",format(CellSize),".tif") #Save the raster writeRaster(NewRaster, filename=OutputFilePath, format="GTiff", overwrite=TRUE) # convert the raster to a matrix that R can perform stats on TheMatrix=as.matrix(NewRaster)
4. Compute the overall mean of the standard deviation raster. This is the average standard deviation for our down-sampled raster.
5. Estimate a mean standard deviation for each pixel in the original raster. A formal process for this would include making measurements in the area the raster originated from and developing an approximation for the standard deviation. We'll estimate this by making a visual examination of the area around campus but if you can visit your study site, you can do this yourself. Remember that the standard deviation contains about 2/3's of the values within one standard deviation of the mean.
6. Write a function in R to inject noise (randomness) into the original raster (you can use the function above as a starting point).
7. Create a loop that injects noise into a raster, finds the overall (mean) standard deviation of the raster, saves this in a vector, and then repeats the process.
8. Plot the histogram of your results.
The image below illustrates this process. The "Original Raster" is the one we started with. We take it and add some "random noise" to create a noisy raster. This is simulating the process of collecting another dataset without having to go out in the field. Typically we would down sample this raster and then use it in our models. However, this exercise is seeing about how noise affects our models so we want to create a raster with the standard deviation of each of the pixels in the Down-sampled raster instead. Then, we'll record the means of the standard deviation rasters and repeat the process a number of times. When done, take a look at the histogram and see what effect the noise is having on your raster covariate layers.
9. Finally, run the loop above with the different cell sizes.
1. Obtain a point data set to create an interpolated surface. If you do not have one, you can use a remotely sensed dataset (such as a DEM) and select points that might be used for interpolation (flat spots, etc.).
2. Create an interpolated surface from the point data set.
3. Find the overall standard deviation for the raster. This will only quantify the uncertainty in the process of interpolating, not the uncertainty from sub-sampling the data or the uncertainty in measurements.
4. Estimate the standard deviation of the data.
5. Randomly subsample the original data set (or raster) and add error based on the standard deviation.
6. Repeat 5 and 6 saving the min, max, and standard deviation values.
Question: What is the effect of different numbers of sub-samples on the uncertainty of our final interpolated raster.
Introduction to the "Raster" package
© Copyright 2018 HSU - All rights reserved.