| Title: | Radiometric and Topographic Correction of Satellite Imagery |
|---|---|
| Description: | Processing of Landsat or other multispectral satellite imagery. Includes relative normalization, image-based radiometric correction, and topographic correction options. The original package description was published as Goslee (2011) <doi:10.18637/jss.v043.i04>, and details of the topographic corrections in Goslee (2012) <doi:10.14358/PERS.78.9.973>. |
| Authors: | Sarah Goslee [aut, cre] |
| Maintainer: | Sarah Goslee <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.1.2 |
| Built: | 2024-11-17 03:09:05 UTC |
| Source: | https://github.com/phiala/landsat |
Finds Bare Soil Line (BSL) and maximum vegetation point.
BSL(band3, band4, method = "quantile", ulimit = 0.99, llimit = 0.005, maxval = 255)BSL(band3, band4, method = "quantile", ulimit = 0.99, llimit = 0.005, maxval = 255)
band3 |
File name or image file (matrix, data frame, or SpatialGridDataFrame) for Landsat band 3 DN (red). |
band4 |
File name or image file (matrix, data frame, or SpatialGridDataFrame) for Landsat band 4 DN (NIR). |
method |
Either "quantile" or "minimum" – describes way in which soil line is identified. |
ulimit |
Upper limit for quantile of band ratios (ulimit < 1). |
llimit |
Lower limit for quantile of band ratios (llimit > 0). |
maxval |
Maximum value for band data; default of 255 for Landsat 5 and 7. |
Finding the BSL requires identifying the lowest NIR values for each level of red. The quantile method takes the lowest set of points, those with a NIR/red ratio less than the llimit-th quantile. The minimum value method takes the lowest NIR value for each level of red. However they are found, these points with low NIR for their red values are used in a major axis regression to find the Bare Soil Line. This function also identifies the full canopy point (maximum vegetation), by using the ulimit to identify the top points, with NIR/red ratio greater than the ulimit-th quantile, and with high NIR values. Red or NIR values of 255 (saturated sensor) are omitted when calculating the BSL.
BSL |
Regression coefficients for the Bare Soil Line |
top |
band 3 and band 4 values for the full canopy point |
Sarah Goslee
Maas, S. J. & Rajan, N. 2010. Normalizing and converting image DC data using scatter plot matching. Remote Sensing 2:1644-1661.
data(nov3) data(nov4) nov.bsl <- BSL(nov3, nov4) plot(as.vector(as.matrix(nov3)), as.vector(as.matrix(nov4))) abline(nov.bsl$BSL, col="red") points(nov.bsl$top[1], nov.bsl$top[2], col="green", cex=2, pch=16)data(nov3) data(nov4) nov.bsl <- BSL(nov3, nov4) plot(as.vector(as.matrix(nov3)), as.vector(as.matrix(nov4))) abline(nov.bsl$BSL, col="red") points(nov.bsl$top[1], nov.bsl$top[2], col="green", cex=2, pch=16)
Uses Landsat band 1 and band 6 to identify clouds and create a cloud mask.
clouds(band1, band6, level = 0.0014, buffer=5)clouds(band1, band6, level = 0.0014, buffer=5)
band1 |
File name or image file (matrix, data frame, or SpatialGridDataFrame) for Landsat band 1. |
band6 |
File name or image file (matrix, data frame, or SpatialGridDataFrame) for Landsat band 6. |
level |
Threshold level for cloud/noncloud decision. The default threshold is appropriate for reflectance and temperature values, and must be adjusted for use with DN. |
buffer |
Pixel buffer size to expand around thresholded cloud areas. |
Clouds are reflective (high) in band 1 and cold (low) in band 6, so the ratio of the two bands is high over clouds. The ratio must be adjusted for data type, whether reflectance, radiance, or DN.
Returns a cloud mask in the same format as band1. Clouds are 1; noncloud areas are NA. Cloud areas are expanded by buffer pixels to ensure that cloud edges are captured.
Sarah Goslee
This function is loosely based on: Martinuzzi, S., Gould, W.A., Ramos Gonzales, O.M. 2007. Creating Cloud-Free Landsat ETM+ Data Sets in Tropical Landscapes: Cloud and Cloud-Shadow Removal. USDA Forest Service General Technical Report IITF-GTR-32.
data(july1) data(july61) july.cloud <- clouds(july1, july61) par(mfrow=c(1,2)) image(july1) image(july.cloud)data(july1) data(july61) july.cloud <- clouds(july1, july61) par(mfrow=c(1,2)) image(july1) image(july.cloud)
Convert a vector containing year, month, day or individual year, month, day arguments into a decimal date in years.
ddate(year, month, day)ddate(year, month, day)
year |
Either a numeric year OR a vector in the form of c(year, month, day). The latter option is so that ddate() can be conveniently used with apply(). |
month |
If year is a single value, must contain the number of the month. |
day |
If year is a single value, must contain the number of the day. |
ddate() will accept a vector with the three date components so that it can be conveniently used with apply() on a data frame containing columns for year, month and day.
The decimal date in years.
Sarah Goslee
ddate(2001, 5, 15)ddate(2001, 5, 15)
A 30-meter resolution elevation model in SpatialGridDataFrame format that matches the Landsat images nov and july.
data(dem)data(dem)
Elevations are in meters.
Digital Elevation Models for the United States are available from the United States Geologic Survey, http://www.usgs.gov
data(dem) dem.slopeasp <- slopeasp(dem) par(mfrow=c(1,3)) image(dem) image(dem.slopeasp$slope) image(dem.slopeasp$aspect)data(dem) dem.slopeasp <- slopeasp(dem) par(mfrow=c(1,3)) image(dem) image(dem.slopeasp$slope) image(dem.slopeasp$aspect)
Calculates calibration value for the Dark Object Subtraction (DOS) method of radiometric correction.
DOS(sat = 5, scattering.coef = c(-4, -2, -1, -0.7, -0.5), SHV, SHV.band, gain, offset, Grescale, Brescale, sunelev, edist, Esun = c(198.3, 179.6, 153.6, 103.1, 22, 8.34), blackadjust = 0.01)DOS(sat = 5, scattering.coef = c(-4, -2, -1, -0.7, -0.5), SHV, SHV.band, gain, offset, Grescale, Brescale, sunelev, edist, Esun = c(198.3, 179.6, 153.6, 103.1, 22, 8.34), blackadjust = 0.01)
sat |
Landsat satellite platform: 5 for TM; 7 for ETM+; 8 for OLI. |
scattering.coef |
Atmospheric scattering coefficient; defaults are from Chavez 1988. |
SHV |
Starting Haze Value |
SHV.band |
Band from which the Starting Haze Value was obtained. |
gain |
Band-specific sensor gain. Require either gain and offset or Grescale and Brescale to convert DN to radiance. |
offset |
Band-specific sensor offset. Require either gain and offset or Grescale and Brescale to convert DN to radiance. |
Grescale |
Band-specific sensor $G_rescale$ (gain). Require either gain and offset or Grescale and Brescale to convert DN to radiance. |
Brescale |
Band-specific sensor $B_rescale$ (bias). Require either gain and offset or Grescale and Brescale to convert DN to radiance. |
sunelev |
Sun elevation in degrees |
edist |
Earth-Sun distance in AU. |
Esun |
Exo-atmospheric solar irradiance, as given by Chander et al. 2009 or others. |
blackadjust |
By default, implements 1% adjustment value to compensate for lack of perfectly dark pixels. |
The Dark Object Subtraction method assumes that the darkest parts of an image (water, artificial structures) should be black if not for the effects of atmospheric scatter. Corrections to make it possible to use the black value from one band to correct the remaining bands.
DNfinal.mean |
The Dark Object Subtraction value for the complete set of scattering coefficients (Table X in Chavez 1989). |
DNfinal.approx |
The same table as DNfinal.mean, but using a numerical approximation across the band wavelength instead of the mean wavelength. |
Sarah Goslee
Chavez, Jr., P. S. 1988. An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data. Remote Sensing of Environment 24:459-479.
Chavez, Jr., P. S. 1989. Radiometric calibration of Landsat Thematic Mapper multispectral images. Photogrammetric Engineering and Remote Sensing 55: 1285-1294.
data(july1) data(july3) # One approach to choosing a Starting Haze Value is to take # the lowest DN value with a frequency greater than some # predetermined threshold, in this case 1000 pixels. SHV <- table(july1@data[,1]) SHV <- min(as.numeric(names(SHV)[SHV > 1000])) # this is used as Lhaze in the radiocorr function # G_rescale, B_rescale, sun elevation comes from metadata for the SHV band july.DOS <- DOS(sat=7, SHV=SHV, SHV.band=1, Grescale=0.77569, Brescale=-6.20000, sunelev=61.4, edist=ESdist("2002-07-20"))$DNfinal.mean # DOS() returns results for the complete set of scattering coefficients # need to choose the appropriate one based on general atmospheric conditions ### -4.0: Very Clear SHV <= 55 ### -2.0: Clear SHV 56-75 ### -1.0: Moderate SHV 76-95 ### -0.7: Hazy SHV 96-115 ### -0.5: Very Hazy SHV >115 # for july, SHV == 70, so use -2.0: Clear july.DOS <- july.DOS[ , 2] # Use DOS value as Lhaze in radiocorr() for DOS correction to reflectance july3.DOSrefl <- radiocorr(july3, Grescale=0.77569, Brescale=-6.20000, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1533, Lhaze=july.DOS[3], method="DOS")data(july1) data(july3) # One approach to choosing a Starting Haze Value is to take # the lowest DN value with a frequency greater than some # predetermined threshold, in this case 1000 pixels. SHV <- table(july1@data[,1]) SHV <- min(as.numeric(names(SHV)[SHV > 1000])) # this is used as Lhaze in the radiocorr function # G_rescale, B_rescale, sun elevation comes from metadata for the SHV band july.DOS <- DOS(sat=7, SHV=SHV, SHV.band=1, Grescale=0.77569, Brescale=-6.20000, sunelev=61.4, edist=ESdist("2002-07-20"))$DNfinal.mean # DOS() returns results for the complete set of scattering coefficients # need to choose the appropriate one based on general atmospheric conditions ### -4.0: Very Clear SHV <= 55 ### -2.0: Clear SHV 56-75 ### -1.0: Moderate SHV 76-95 ### -0.7: Hazy SHV 96-115 ### -0.5: Very Hazy SHV >115 # for july, SHV == 70, so use -2.0: Clear july.DOS <- july.DOS[ , 2] # Use DOS value as Lhaze in radiocorr() for DOS correction to reflectance july3.DOSrefl <- radiocorr(july3, Grescale=0.77569, Brescale=-6.20000, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1533, Lhaze=july.DOS[3], method="DOS")
Calculates the estimated Earth-Sun distance in Astronomical Units (AU) for a given date.
ESdist(adate)ESdist(adate)
adate |
date in "YYYY-MM-DD" format |
Returns estimated Earth-Sun distance in AU.
Sarah Goslee
ESdist("2010-08-30")ESdist("2010-08-30")
Finds best fit between target image and tofix image by minimizing RMSE between the two. The tofix image is moved one pixel at a time horizontally or vertically. Simple automated georeferencing is adequate for some image-processing tasks.
georef(target, tofix, maxdist = 1000, startx = 0, starty = 0)georef(target, tofix, maxdist = 1000, startx = 0, starty = 0)
target |
A georeferenced base image; can be matrix, dataframe or SpatialGridDataFrame. |
tofix |
The image to be georeferenced; can be matrix, dataframe or SpatialGridDataFrame. |
maxdist |
The greatest distance to move the tofix image. If this is exceeded, the function will stop. |
startx |
Shift the tofix image this many pixels in the x direction before beginning, to avoid local minimum. |
starty |
Shift the tofix image this many pixels in the y direction before beginning, to avoid local minimum. |
This function offers a simplistic approach to georeferencing using an iterative algorithm that at each step moves the tofix image one pixel in the direction that produces the greatest reduction in RMSE. When RMSE no longer decreases or maxdist is reached, the algorithm stops, assuming that the tofix image now matches the reference target image. This algorithm can produce local minima. Results should always be checked visually.
Note: this algorithm is only effective with images larger than the samples included with this package. The July and November images are already georectified, but this function will show them as needing considerable adjustment. Images of at least 1000x1000 pixels are necessary for adequate results.
shiftx |
The x-direction shift to get the best match (lowest RMSE). |
shifty |
The y-direction shift to get the best match (lowest RMSE). |
initrmse |
Initial RMSE between target and tofix images. |
currrmse |
Lowest RMSE, after shiftx and shifty pixel adjustments. Will be 9999 if maxdist is exceeded. |
Sarah Goslee
# to use for georeferencing data(nov3) data(july3) july.shift <- georef(nov3, july3, maxdist=50) # match july to november july3.corr <- geoshift(july3, padx=50, pady=50, july.shift$shiftx, july.shift$shifty) # only need to run georef once for a particular date # use the same correction for all bands data(july4) july4.corr <- geoshift(july4, padx=50, pady=50, july.shift$shiftx, july.shift$shifty)# to use for georeferencing data(nov3) data(july3) july.shift <- georef(nov3, july3, maxdist=50) # match july to november july3.corr <- geoshift(july3, padx=50, pady=50, july.shift$shiftx, july.shift$shifty) # only need to run georef once for a particular date # use the same correction for all bands data(july4) july4.corr <- geoshift(july4, padx=50, pady=50, july.shift$shiftx, july.shift$shifty)
Shifts an image vertically or horizontally and adds a padded border.
geoshift(mat, padx, pady, shiftx, shifty, nodata = NA)geoshift(mat, padx, pady, shiftx, shifty, nodata = NA)
mat |
A matrix, data frame or SpatialGridDataFrame |
padx |
Number of pixels to add as padding in the x direction on each side of the image (along the x-axis). Should be larger than the number of pixels to shift to avoid data loss. |
pady |
Number of pixels to add as padding in the y direction on each side of the image (along the y-axis). Should be larger than the number of pixels to shift to avoid data loss. |
shiftx |
Number of pixels to shift (positive or negative) in the x direction (along the x-axis). |
shifty |
Number of pixels to shift (positive or negative) in the y direction (along the y-axis). |
nodata |
Value to use for missing data. |
This function can be used to correct spatially-referenced images that are off by a few pixels in the x or y directions. It does not warp an image, only slide it. Adding padding to the outside edge makes it possible to match several images even if they are not stored with georeferecing information. geoshift() can be used in conjunction with georef() to automatically match up geospatial images. Note: directions are relative to the image as displayed by the image() command, and not the underlying matrix representation.
Returns data in the same format as the function was given: matrix, data frame, or SpatialGridDataFrame.
Sarah Goslee
testmat <- matrix(1:9, 3, 3) geoshift(testmat, 5, 10, 0, 0) geoshift(testmat, 5, 10, 2, 2) # to use for georeferencing data(nov3) data(july3) july.shift <- georef(nov3, july3, maxdist=50) # match july to november july3.corr <- geoshift(july3, padx=50, pady=50, july.shift$shiftx, july.shift$shifty) # only need to run georef once for a particular date # use the same correction for all bands data(july4) july4.corr <- geoshift(july4, padx=50, pady=50, july.shift$shiftx, july.shift$shifty)testmat <- matrix(1:9, 3, 3) geoshift(testmat, 5, 10, 0, 0) geoshift(testmat, 5, 10, 2, 2) # to use for georeferencing data(nov3) data(july3) july.shift <- georef(nov3, july3, maxdist=50) # match july to november july3.corr <- geoshift(july3, padx=50, pady=50, july.shift$shiftx, july.shift$shifty) # only need to run georef once for a particular date # use the same correction for all bands data(july4) july4.corr <- geoshift(july4, padx=50, pady=50, july.shift$shiftx, july.shift$shifty)
Force image x to match target image by matching their histograms.
histmatch(master, tofix, mask, minval = 0, maxval = 255, by = 1)histmatch(master, tofix, mask, minval = 0, maxval = 255, by = 1)
master |
The target image, in SpatialGridDataFrame, data frame, matrix or vector format. |
tofix |
The image to be normalized, in any format. |
mask |
Areas to be omitted, if any, such as a cloud mask. Only NA values within the mask will be used. |
minval |
Lower bound of the possible range of values in target and tofix images. |
maxval |
Upper bound of the possible range of values in target and tofix images. |
by |
Step size to use in constructing histograms. Should be appropriate for minval and maxval of the images. |
The histogram of the tofix image will be forced to match that of the target image.
recode |
The transformation table used to match the histograms. |
newimage |
The transformed image, in the same format in which tofix was provided. |
Sarah Goslee
## Not run: data(nov3) data(july3) par(mfrow=c(2,2)) image(nov3) image(july3) nov3.newR <- relnorm(master=july3, tofix=nov3) image(nov3.newR$newimage) nov3.newH <- histmatch(master=july3, tofix=nov3) image(nov3.newH$newimage) ## End(Not run)## Not run: data(nov3) data(july3) par(mfrow=c(2,2)) image(nov3) image(july3) nov3.newR <- relnorm(master=july3, tofix=nov3) image(nov3.newR$newimage) nov3.newH <- histmatch(master=july3, tofix=nov3) image(nov3.newH$newimage) ## End(Not run)
SpatialGridDataFrame containing a 300 x 300 pixel subset (1500 x 1500 m) of the Landsat ETM+ image for path 15, row 32, obtained on 20 July 2002. Each band, including both thermal bands, is contained in a separate file.
data(july1)data(july1)
Images are in SpatialGridDataFrame format. More information is available in the documentation for the sp package.
Date: 2002-07-20
Satellite: Landsat ETM+ (7)
Sun elevation: 61.4
Sun azimuth: 125.8
——–
band Grescale Brescale
1 0.77569 -6.20
2 0.79569 -6.40
3 0.61922 -5.00
4 0.63725 -5.10
5 0.12573 -1.00
7 0.04373 -0.35
Landsat images can be obtained from the United States Geological Survey at http://landsat.usgs.gov
data(july3) image(july3)data(july3) image(july3)
Uses GDAL tools to reproject (optional) and subset a geotiff given the center point and the desired size.
lssub(filename, outname, centerx, centery, centerepsg, widthx, widthy)lssub(filename, outname, centerx, centery, centerepsg, widthx, widthy)
filename |
Filename (and path) to a geotiff image. |
outname |
Filename (and path) for subset image. |
centerx |
x coordinate of new center point. |
centery |
y coordinate of new center point. |
centerepsg |
Projection of the center point coordinates as 5-digit EPSG code. If missing, assume that point and geotiff have the same projection. |
widthx |
Desired width of subset image. |
widthy |
Desired height of subset image. |
The new image will be a subset of size (widthx, widthy) with center point (centerx, centery), with the same pixel size. If the center point coordinates are in a different projection than the original image, they will be reprojected.
The new image is exported as a geotiff. Nothing is returned within R.
Requires gdalinfo and gdaltransform to be available to the operating system. Only known to work on linux. This function was written to speed processing of multiple files for a specific project, and may be dropped in future releases of the landsat package. On my computer, lssub() is over an order of magnitude faster than reading the image into R, subsetting it, and writing out the result.
Sarah Goslee
## Not run: lssub("/data/gis/testimage.tif", "/data/gis/subimage.tif", centerx = 260485, centery = 4527220, centerepsg = 26918, widthx = 50, widthy = 50) ## End(Not run)## Not run: lssub("/data/gis/testimage.tif", "/data/gis/subimage.tif", centerx = 260485, centery = 4527220, centerepsg = 26918, widthx = 50, widthy = 50) ## End(Not run)
Adds several modified Minnaert corrections to the capabilities of topocorr().
minnaert(x, slope, aspect, sunelev, sunazimuth, na.value = NA, GRASS.aspect=FALSE, IL.epsilon=0.000001, slopeclass = c(1, 5, 10, 15, 20, 25, 30, 45), coverclass)minnaert(x, slope, aspect, sunelev, sunazimuth, na.value = NA, GRASS.aspect=FALSE, IL.epsilon=0.000001, slopeclass = c(1, 5, 10, 15, 20, 25, 30, 45), coverclass)
x |
Image to be corrected, in matrix, data frame, or SpatialGridDataFrame format. |
slope |
Slope image of same size and resolution as x. |
aspect |
Aspect image of same size and resolution as x. |
sunelev |
Sun elevation in degrees. |
sunazimuth |
Sun azimuth in degrees. |
na.value |
Value to use for missing data. |
GRASS.aspect |
Whether aspect is measured according to GRASS defaults (counterclockwise from east) or is measured clockwise from north. If GRASS.aspect=TRUE, aspect is converted to clockwise from north before analysis. |
IL.epsilon |
If IL == 0 (Illumination), some methods will give a topographically-corrected value of Inf due to division by zero. If desired, adding a small increment to zero values eliminates this. |
slopeclass |
The classes into which the slope will be divided before calculating k separately for each class. |
coverclass |
If present, TRUE/FALSE vector indicating which pixels to use when calculating k. This allows k to be determined separately for different cover classes. |
Calculates the Minnaert k coefficients for the whole image and for the individual slope classes.
allcoef |
The Minnaert k for the entire image. This is the value used in topocorr() (though the latter may have been truncated). |
classcoef |
A data frame containing the slope class midpoints, number of pixels per class, and k for that class (for the desired cover class, if specified). |
xout |
A topographically-corrected image in the same format as x. |
xout |
A topographically-corrected image in the same format as x. |
Sarah Goslee
Lu, D., Ge, H., He, S., Xu, A., Zhou, G., and Du, H. 2008. Pixel-based Minnaert correction method for reducing topographic effects on a Landsat 7 ETM+ image. Photogrammetric Engineering and Remote Sensing 74:1343-1350.
# require slope and aspect for topographic correction data(dem) dem.slopeasp <- slopeasp(dem) # use cosine method of topographic correction data(july4) july4.minpix <- minnaert(july4, dem.slopeasp$slope, dem.slopeasp$aspect, sunelev=61.4, sunazimuth=125.8, slopeclass=c(1, 5, 10, 15, 50)) july4.minpix$classcoef # all coefficients# require slope and aspect for topographic correction data(dem) dem.slopeasp <- slopeasp(dem) # use cosine method of topographic correction data(july4) july4.minpix <- minnaert(july4, dem.slopeasp$slope, dem.slopeasp$aspect, sunelev=61.4, sunazimuth=125.8, slopeclass=c(1, 5, 10, 15, 50)) july4.minpix$classcoef # all coefficients
Very simple function to apply a kernel to a matrix across a moving window.
movingwindow(x, kernel, na.rm=TRUE)movingwindow(x, kernel, na.rm=TRUE)
x |
A matrix. |
kernel |
The kernel to be applied to the matrix, for example a Sobel kernel. |
na.rm |
NA handling option to be passed to sum(). If TRUE, NA will be returned if any value under the kernel is NA or NaN, otherwise NA values will be omitted. |
This function is used in the calculation of slope and aspect by slopeasp().
Returns the transformed matrix.
Should be rewritten in C for greater efficiency.
Sarah Goslee
data(dem) dem.smoothed <- movingwindow(dem, matrix(c(1,1,1,1,0,1,1,1,1), 3, 3)/8) par(mfrow=c(1,2)) image(dem) image(dem.smoothed)data(dem) dem.smoothed <- movingwindow(dem, matrix(c(1,1,1,1,0,1,1,1,1), 3, 3)/8) par(mfrow=c(1,2)) image(dem) image(dem.smoothed)
SpatialGridDataFrame containing a 300 x 300 pixel subset (1500 x 1500 m) of the Landsat ETM+ image for path 15, row 32, obtained on 25 November 2002. Each band, including both thermal bands, is contained in a separate file.
data(nov1)data(nov1)
Images are in SpatialGridDataFrame format. More information is available in the documentation for the sp package.
Date: 2002-11-25
Satellite: Landsat ETM+ (7)
Sun elevation: 26.2
Sun azimuth: 159.5
——–
band Grescale Brescale
1 0.77569 -6.20
2 0.79569 -6.40
3 0.61922 -5.00
4 0.63725 -5.10
5 0.12573 -1.00
7 0.04373 -0.35
Landsat images can be obtained from the United States Geological Survey at http://landsat.usgs.gov
data(nov3) image(nov3)data(nov3) image(nov3)
Pseudo-invariant features identification for relative radiometric normalization.
PIF(band3, band4, band7, level = 0.99)PIF(band3, band4, band7, level = 0.99)
band3 |
Landsat band 3, as a filename to be imported, a matrix, data frame, or SpatialGridDataFrame. |
band4 |
Landsat band 4, as a filename to be imported, a matrix, data frame, or SpatialGridDataFrame. |
band7 |
Landsat band 7, as a filename to be imported, a matrix, data frame, or SpatialGridDataFrame. |
level |
Threshold level for identifying PIFs. (0 < level < 1) |
Pseudo-invariant features (PIFs) are areas such as artificial structures that can reasonably be expected to have a constant reflectance over time, rather than varying seasonally as vegetation does. Differences in PIF reflectance between dates can be assumed to be due to varying atmospheric conditions.
Returns a PIF mask in the same format as the input files, with 1 for pseudo-invariant features and 0 for background data.
Sarah Goslee
Schott, J. R.; Salvaggio, C. & Volchok, W. J. 1988. Radiometric scene normalization using pseudoinvariant features. Remote Sensing of Environment 26:1-16.
# identify pseudo-invariant feature data(july3) data(july4) data(july7) july.pif <- PIF(july3, july4, july7) # use PIFs to related nov to july Landsat data for band 3 # properly, would also remove cloudy areas first data(nov3) # use major axis regression: error in both x and y nov.correction <- lmodel2:::lmodel2(july3@data[july.pif@data[,1] == 1, 1] ~ nov3@data[july.pif@data[,1] == 1, 1])$regression.results[2, 2:3] nov3.corrected <- nov3 nov3.corrected@data[,1] <- nov3@data[,1] * nov.correction[2] + nov.correction[1]# identify pseudo-invariant feature data(july3) data(july4) data(july7) july.pif <- PIF(july3, july4, july7) # use PIFs to related nov to july Landsat data for band 3 # properly, would also remove cloudy areas first data(nov3) # use major axis regression: error in both x and y nov.correction <- lmodel2:::lmodel2(july3@data[july.pif@data[,1] == 1, 1] ~ nov3@data[july.pif@data[,1] == 1, 1])$regression.results[2, 2:3] nov3.corrected <- nov3 nov3.corrected@data[,1] <- nov3@data[,1] * nov.correction[2] + nov.correction[1]
Implements several different methods for absolute radiometric correction of satellite data.
radiocorr(x, gain, offset, Grescale, Brescale, sunelev, satzenith = 0, edist, Esun, Lhaze, method = "apparentreflectance")radiocorr(x, gain, offset, Grescale, Brescale, sunelev, satzenith = 0, edist, Esun, Lhaze, method = "apparentreflectance")
x |
Image to be corrected, in matrix, data frame, or SpatialGridDataFrame format. |
gain |
Band-specific sensor gain. Require either gain and offset or Grescale and Brescale to convert DN to radiance. |
offset |
Band-specific sensor offset. Require either gain and offset or Grescale and Brescale to convert DN to radiance. |
Grescale |
Band-specific sensor Grescale (gain). Require either gain and offset or Grescale and Brescale to convert DN to radiance. |
Brescale |
Band-specific sensor Brescale (bias). Require either gain and offset or Grescale and Brescale to convert DN to radiance. |
sunelev |
Sun elevation in degrees |
satzenith |
Satellite sensor zenith angle (0 for Landsat) |
edist |
Earth-Sun distance in AU. |
Esun |
Exo-atmospheric solar irradiance, as given by Chandler et al. 2009 or others. |
Lhaze |
Haze value, such as SHV from DOS() function. Not needed for apparent reflectance. |
method |
Radiometric correction method to be used. There are currently four methods available: "apparentreflectance", "DOS" (Chavez 1989), "COSTZ" (Chavez 1996), "DOS4" (SWS+2001). |
Uses one of four image-based radiometric correction methods to adjust a satellite image to compensate for atmospheric conditions.
Returns a radiometrically-corrected image in the same format as x.
Sarah Goslee
Chavez, Jr., P. S. 1989. Radiometric calibration of Landsat Thematic Mapper multispectral images. Photogrammetric Engineering and Remote Sensing 55:1285-1294.
Chavez, Jr., P. S. 1996. Image-based atmospheric corrections revisited and improved. Photogrammetric Engineering and Remote Sensing 62:1025-1036.
Song, C.; Woodcock, C. E.; Seto, K. C.; Lenney, M. P. & Macomber, S. A. 2001. Classification and change detection using Landsat TM data: when and how to correct atmospheric effects? Remote Sensing of Environment 75:230-244.
data(july1) data(july3) # One approach to choosing a Starting Haze Value is to take the lowest DN value # with a frequency greater than some predetermined threshold, in this case 1000 pixels. SHV <- table(july1@data[,1]) SHV <- min(as.numeric(names(SHV)[SHV > 1000])) # this is used as Lhaze in the radiocorr function # Grescale, Brescale, sun elevation comes from metadata for the SHV band july.DOS <- DOS(sat=7, SHV=SHV, SHV.band=1, Grescale=0.77569, Brescale=-6.20000, sunelev=61.4, edist=ESdist("2002-07-20"))$DNfinal.mean # DOS() returns results for the complete set of scattering coefficients # need to choose the appropriate one based on general atmospheric conditions ### -4.0: Very Clear SHV <= 55 ### -2.0: Clear SHV 56-75 ### -1.0: Moderate SHV 76-95 ### -0.7: Hazy SHV 96-115 ### -0.5: Very Hazy SHV >115 # for july, SHV == 70, so use -2.0: Clear july.DOS <- july.DOS[ , 2] # Use DOS value as Lhaze in radiocorr() for DOS correction to reflectance july3.DOSrefl <- radiocorr(july3, Grescale=0.77569, Brescale=-6.20000, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1533, Lhaze=july.DOS[3], method="DOS")data(july1) data(july3) # One approach to choosing a Starting Haze Value is to take the lowest DN value # with a frequency greater than some predetermined threshold, in this case 1000 pixels. SHV <- table(july1@data[,1]) SHV <- min(as.numeric(names(SHV)[SHV > 1000])) # this is used as Lhaze in the radiocorr function # Grescale, Brescale, sun elevation comes from metadata for the SHV band july.DOS <- DOS(sat=7, SHV=SHV, SHV.band=1, Grescale=0.77569, Brescale=-6.20000, sunelev=61.4, edist=ESdist("2002-07-20"))$DNfinal.mean # DOS() returns results for the complete set of scattering coefficients # need to choose the appropriate one based on general atmospheric conditions ### -4.0: Very Clear SHV <= 55 ### -2.0: Clear SHV 56-75 ### -1.0: Moderate SHV 76-95 ### -0.7: Hazy SHV 96-115 ### -0.5: Very Hazy SHV >115 # for july, SHV == 70, so use -2.0: Clear july.DOS <- july.DOS[ , 2] # Use DOS value as Lhaze in radiocorr() for DOS correction to reflectance july3.DOSrefl <- radiocorr(july3, Grescale=0.77569, Brescale=-6.20000, sunelev=61.4, edist=ESdist("2002-07-20"), Esun=1533, Lhaze=july.DOS[3], method="DOS")
The Radiometric Control Sets method of relative radiometric correction for Landsat data.
RCS(data.tc, level = 0.01)RCS(data.tc, level = 0.01)
data.tc |
The output of tasscap(). |
level |
Threshold level to use (0 < level < 1). |
Radiometric Control Sets (RCSs) are areas such as artificial structures and large bodies of water that can reasonably be expected to have a constant reflectance over time, rather than varying seasonally as vegetation does. Differences in RCS reflectance between dates can be assumed to be due to varying atmospheric conditions. Pixels with low greenness and either high or low brightness are identified.
Returns an RCS mask file in the format of the original data (vector, matrix, data frame or SpatialGridDataFrame, as preseved by tasscap()) with 1 for RCS pixels and 0 for background.
Sarah Goslee
Hall, F.; Strebel, D.; Nickeson, J. & Goetz, S. 1991. Radiometric rectification: toward a common radiometric response among multidate, multisensor images. Remote Sensing of Environment 35:11-27.
# identify radiometric control set data(july1) data(july2) data(july3) data(july4) data(july5) data(july7) july.tc <- tasscap("july", 7) july.rcs <- RCS(july.tc) # use RCS to relate nov to july Landsat data for band 3 # properly, would also remove cloudy areas first data(nov3) # use major axis regression: error in both x and y nov.correction <- lmodel2:::lmodel2(july3@data[july.rcs@data[,1] == 1, 1] ~ nov3@data[july.rcs@data[,1] == 1, 1])$regression.results[2, 2:3] nov3.corrected <- nov3 nov3.corrected@data[,1] <- nov3@data[,1] * nov.correction[2] + nov.correction[1]# identify radiometric control set data(july1) data(july2) data(july3) data(july4) data(july5) data(july7) july.tc <- tasscap("july", 7) july.rcs <- RCS(july.tc) # use RCS to relate nov to july Landsat data for band 3 # properly, would also remove cloudy areas first data(nov3) # use major axis regression: error in both x and y nov.correction <- lmodel2:::lmodel2(july3@data[july.rcs@data[,1] == 1, 1] ~ nov3@data[july.rcs@data[,1] == 1, 1])$regression.results[2, 2:3] nov3.corrected <- nov3 nov3.corrected@data[,1] <- nov3@data[,1] * nov.correction[2] + nov.correction[1]
Use regression methods to adjust distribution of values in image tofix to match those in the master image.
relnorm(master, tofix, mask, method = "MA", nperm = 1000)relnorm(master, tofix, mask, method = "MA", nperm = 1000)
master |
The target image, in SpatialGridDataFrame, data frame, matrix or vector format. |
tofix |
The image to be normalized, in any format. |
mask |
Areas to be omitted, if any, such as a cloud mask. Only NA values within the mask will be used. |
method |
Regression method to be used. OLS: Ordinary Least Squares; MA: Major Axis (recommended); SMA: Standard Major Axis. |
nperm |
Number of permutations to use for significance testing. |
The regression coefficients from tofix ~ master will be used to match the distribution of values of tofix to those in the master image.
regression.results |
The regression results from lmodel2 |
newimage |
The transformed image, in the same format in which tofix was provided. |
Sarah Goslee
## Not run: data(nov3) data(july3) par(mfrow=c(2,2)) image(nov3) image(july3) nov3.newR <- relnorm(master=july3, tofix=nov3) image(nov3.newR$newimage) nov3.newH <- histmatch(master=july3, tofix=nov3) image(nov3.newH$newimage) ## End(Not run)## Not run: data(nov3) data(july3) par(mfrow=c(2,2)) image(nov3) image(july3) nov3.newR <- relnorm(master=july3, tofix=nov3) image(nov3.newR$newimage) nov3.newH <- histmatch(master=july3, tofix=nov3) image(nov3.newH$newimage) ## End(Not run)
Uses gridded elevation data to calculate slope and aspect, by default using a 3x3 region. The horizontal resolution and vertical resolution must be in the same units.
slopeasp(x, EWres, NSres, EWkernel, NSkernel, smoothing=1)slopeasp(x, EWres, NSres, EWkernel, NSkernel, smoothing=1)
x |
gridded elevation data, either as a SpatialGridDataFrame, dataframe, or matrix. |
EWres |
East-West grid resolution. May be omitted if x is a SpatialGridDataFrame and the horizontal units are the same as the vertical units. |
NSres |
North-South grid resolution. May be omitted if x is a SpatialGridDataFrame and the horizontal units are the same as the vertical units. |
EWkernel |
The kernel to use when calculating the East-West component of slope. If missing, a 3x3 kernel will be used. |
NSkernel |
The kernel to use when calculating the North-South component of slope. If missing, a 3x3 kernel will be used. |
smoothing |
A positive integer describing the additional smoothing to be applied, if any. smoothing=1 (default) means no smoothing will be used. |
By default, a 3x3 Sobel filter is used (as is standard in many GIS packages). A larger Sobel filter or a different filter will give varying results. This filter provides the third-order finite difference weighted by reciprocal of distance method proposed by Unwin (1981).
slope |
The slope of the DEM, in degrees |
aspect |
The aspect of the DEM, beginning with north and moving clockwise, and with aspect = 0 where slope = 0. |
Sarah Goslee
Unwin, D. 1981. Introductory Spatial Analysis. London: Methuen. Clarke, K.C.and Lee, S.J. 2007. Spatial resolution and algorithm choice as modifiers of downslope flow computed from Digital Elevation Models. Cartography and Geographic Information Science 34:215-230.
data(dem) dem.slopeasp <- slopeasp(dem) par(mfrow=c(1,3)) image(dem) image(dem.slopeasp$slope) image(dem.slopeasp$aspect)data(dem) dem.slopeasp <- slopeasp(dem) par(mfrow=c(1,3)) image(dem) image(dem.slopeasp$slope) image(dem.slopeasp$aspect)
Tasseled cap transformation for Landsat TM, ETM+, or OLI.
tasscap(basename, sat = 7)tasscap(basename, sat = 7)
basename |
Base filename (string) to which band numbers are appended, eg "july" for files named "july1", "july2", "july3", etc. Data should be at-sensor reflectance. |
sat |
Landsat satellite platform: 5 for TM; 7 for ETM+; 8 for OLI. |
For Landsat TM, the coefficients are to be applied to "reflectance factors", which appear to be the DN. For ETM+ and OLI, the coefficients are for top-of-atmosphere reflectance. For both TM and ETM+, the bands to be provided are 1, 2, 3, 4, 5, and 7. For OLI, the bands needed are 2 through 7. Future updates will allow use of a raster stack rather than separate objects.
If the input files are matrices or data frames, returns a data frame with three columns, one for each component. If the input files are SpatialGridDataFrames, returns a list with one element for each component. In either case three components are returned: Brightness, Greenness, Wetness.
Sarah Goslee
Original papers:
Baig, M. H. A., Zhang, L., Shuai, T. & Tong, Q. 2014. Derivation of a tasselled cap transformation based on Landsat 8 at-satellite reflectance. Remote Sensing Letters 5:423-431.
Crist, E. P. 1985. A TM tasseled cap equivalent transformation for reflectance factor data. Remote Sensing of Environment 17:301-306.
Crist, E. & Kauth, R. 1986. The tasseled cap de-mystified. Photogrammetric Engineering and Remote Sensing 52:81-86.
Huang, C., Wylie, B., Yang, L., Homer, C. & Zylstra, G. 2002. Derivation of a tasseled cap transformation based on Landsat 7 at-satellite reflectance. International Journal of Remote Sensing 23:1741-1748.
data(july1) data(july2) data(july3) data(july4) data(july5) data(july7) july.tc <- tasscap("july", 7)data(july1) data(july2) data(july3) data(july4) data(july5) data(july7) july.tc <- tasscap("july", 7)
Converts Landsat thermal band DN (TM or ETM+ band 6-1 and 6-2) to temperature using default coefficients from Chander et al. 2009.
thermalband(x, band)thermalband(x, band)
x |
Landsat band 6 Digital Number (DN) in matrix, data frame or SpatialGridDataFrame format. |
band |
6 for TM; 61 or 62 for the appropriate ETM+ bands. Any other value will fail. |
Returns a temperature image in the same format as x.
Sarah Goslee
Coefficients from Chander, G., Markham, B.L., Helder, D.L. 2009. Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote Sensing of Environment 113:893-903.
data(nov61) nov.temp1 <- thermalband(nov61, 61) image(nov.temp1)data(nov61) nov.temp1 <- thermalband(nov61, 61) image(nov.temp1)
Implements several different methods for topographic correction of remote sensing data.
topocorr(x, slope, aspect, sunelev, sunazimuth, method = "cosine", na.value = NA, GRASS.aspect=FALSE, IL.epsilon=0.000001)topocorr(x, slope, aspect, sunelev, sunazimuth, method = "cosine", na.value = NA, GRASS.aspect=FALSE, IL.epsilon=0.000001)
x |
Image to be corrected, in matrix, data frame, or SpatialGridDataFrame format. |
slope |
Slope image of same size and resolution as x. |
aspect |
Aspect image of same size and resolution as x. |
sunelev |
Sun elevation in degrees. |
sunazimuth |
Sun azimuth in degrees. |
method |
Topographic correction method to be used. There are currently eight methods available: "cosine", "improvedcosine", "minnaert", "ccorrection" (first four from Riano et al. 2003), "minslope" (Minnaert with slope correction, also from Riano et al. 2003), "gamma" (from Richter et al. 2009), "SCS" (Gu and Gillespie 1998, Gao and Zhang 2009), "illumination" (uncorrected illumination). |
na.value |
Value to use for missing data. |
GRASS.aspect |
Whether aspect is measured according to GRASS defaults (counterclockwise from east) or is measured clockwise from north. If GRASS.aspect=TRUE, aspect is converted to clockwise from north before analysis. |
IL.epsilon |
If IL == 0 (Illumination), some methods will give a topographically-corrected value of Inf due to division by zero. If desired, adding a small increment to zero values eliminates this. |
Uses one of the available topographic correction methods to compensate for the effects of slope and aspect on reflectance from the land surface.
Returns a topographically-corrected image in the same format as x.
Sarah Goslee
Gao, Y. & Zhang, W. 2009. LULC classification and topographic correction of Landsat-7 ETM+ imagery in the Yangjia River Watershed: the influence of DEM resolution. Sensors 9:1980-1995.
Gu, D. & Gillespie, A. 1998. Topographic normalization of Landsat TM images of forest based on subpixel sun-canopy-sensor geometry. Remote Sensing of Environment 64:166-175.
Riano, D., Chuvieco, E., Salas, J. & Aguado, I. 2003. Assessment of different topographic corrections in Landsat-TM data for mapping vegetation types. IEEE Transactions on Geoscience and Remote Sensing 41:1056-1061.
Richter, R., Kellenberger, T. & Kaufmann, H. 2009. Comparison of topographic correction methods. Remote Sensing 1:184-196.
# require slope and aspect for topographic correction data(dem) dem.slopeasp <- slopeasp(dem) # use cosine method of topographic correction data(july3) july3.topo <- topocorr(july3, dem.slopeasp$slope, dem.slopeasp$aspect, sunelev=61.4, sunazimuth=125.8)# require slope and aspect for topographic correction data(dem) dem.slopeasp <- slopeasp(dem) # use cosine method of topographic correction data(july3) july3.topo <- topocorr(july3, dem.slopeasp$slope, dem.slopeasp$aspect, sunelev=61.4, sunazimuth=125.8)