Compare two rasters with raster overlay and replace cell value with variable

I want to compare 2 rasters and when there is a change from 2 to 1 create a new raster that replaces the cell where the change has been detected with the change variable.

## Finds Cells that have changes from 2 to 1 Test <- (rasterT1==2 & rasterT2==1) Date <- Test Date[Test == 1,] <- change

Essentially is if rasterT1 > rasterT2 then add the change value to that cell. With the new raster called year using raster overlay, however I'm not sure how this would be executed.

If you need any more details let me know.

A simple "ifelse" statement should suffice in evaluating a condition.

Here we create two random vectors of [1,2] and apply an ifelse to evaluate the condition of if x = 2 and y = 1 THEN change (1) ELSE no change (0).

( x <- round(runif(10, 1, 2)) ) ( y <- round(runif(10, 1, 2)) ) ifelse( x == 2 & y == 1, 1, 0)

Since this is just an evaluation of one vector being greater than another then one can just use an operator and coerce the Boolean result to numeric [0,1] corresponding to FALSE/TRUE.

as.numeric(x > y)

We can then expand this application of ifelse to a raster stack using overlay. It is best to encode change with an actual value rather than a character. The resulting raster will have the values, as in the above example, 1 for change and 0 for no change.

library(raster) fn <- system.file("external/test.grd", package="raster") s <- stack(fn, fn) s[[1]] <- round(runif(ncell(s), 1, 2)) s[[2]] <- round(runif(ncell(s), 1, 2)) s.change <- overlay(s[[1]], s[[2]], fun = function(x,y) { ifelse( x == 2 & y == 1, 1, 0) } )

If it is always going to be a "> or <" operator, then things are much simpler. The double brackets index a specific raster in the stack.

s.change <- s[[1]] > s[[2]]

Bring rasters to same extent by filling with NA - in R

I have several cropped rasters with different geometry/outlines. Specifically spatial yield maps from several years of the same field, but the extent varies - the measurements were not always overall the whole field, but in some years only part of it. I want to calculate a mean value of those maps and combine them into one mean-value-raster. That does mean however, that not for every pixel in let's say 5 layers/rasters there is a value. I could accept these missing values to be NA, so the final mean value would only be calculated by let's say 3 rasters for parts of the field, where the maps are not overlapping.

I thought of extending the raster with 'extend', filling the non-overlapping parts with NA values:

y <- extend(y, shape, value=NA) #Shape is a rectangular shape that enframes all yield map rasters

That works fine, for all rasters. But they still don't have the same extent. Even if I adjust the extent by setExtent() or extent() <- extent() to the extent of the rectangular shapefile or even to one of the other extended rasters, I still get:

Error in compareRaster(x) : different number or columns

..when I want to stack them and use calc(y, fun=mean. ) . The original raster extents are too different to resample. But they do have the same resolution and CRS.

You can just use min with two raster layers.
Let's start with a reproducible example:

Now we calculate the min of each overlapping cell within the two raster layers very easily with the implemented cell statistics:

Furthermore, you can also apply statistics like mean , max , etc.

If the implemented statistics somehow fail, or you want to use your own statistics, you can also directly access the data per pixel. That is, you first copy one of the raster layers.

Calc : Calculate

Calculate values for a new Raster* object from another Raster* object, using a formula.

If x is a RasterLayer, fun is typically a function that can take a single vector as input, and return a vector of values of the same length (e.g. sqrt ). If x is a RasterStack or RasterBrick, fun should operate on a vector of values (one vector for each cell). calc returns a RasterLayer if fun returns a single value (e.g. sum ) and it returns a RasterBrick if fun returns more than one number, e.g., fun=quantile .

In many cases, what can be achieved with calc , can also be accomplished with a more intuitive 'raster-algebra' notation (see Arith-methods ). For example, r <- r * 2 instead of

r <- calc(r, fun=function(x) , or r <- sum(s) instead of

r <- calc(s, fun=sum) . However, calc should be faster when using complex formulas on large datasets. With calc it is possible to set an output filename and file type preferences.

See ( overlay ) to use functions that refer to specific layers, like ( function(a,b,c) )

Remote Sensing and GIScience in Geomorphology

M.J. Smith , . M. Geilhausen , in Treatise on Geomorphology , 2013 Display

Raster outputs of terrain are usually displayed and inspected on a video display unit (VDU) using the red, green, and blue (RGB) additive model of mixing colors ( Figure 12 ). Given the sensitivity of the human eye, the color cube allows the mixing of the three primary colors at different intensities to provide the full gamut of possible viewable colors. It is therefore possible to display up to three separate images of the same area at any one time, by encoding the pixel values as a 3D coordinate position within the color cube. This is known as a ‘false-color composite’ (unless it represents the actual object color as perceived by the human eye in which case it is a ‘true-color composite’). The technique is widely used in remote sensing and, within geovisualization, allows the interactive inspection of multiple terrain datasets.

Figure 12 . A red, green, and blue (RGB) color cube. The cube defines the intensity of each primary color, allowing the 3D plotting of all colors within the cube. The gray line shows how colors are plotted for a gray scale each where each RGB component has equal intensity.

Figure 13 presents an example of how a false-color composite can be used to interpret the landscape, as it combines slope angle (red: 3×3 kernel), topographic openness (green: 51×51 kernel), and topographic openness (blue: 501×501 kernel). Reddish colors represent high slope and low openness (e.g., at steep and enclosed hillslopes), whereas greenish/bluish colors depict low slopes and relatively high openness (e.g., in flat and open areas).

Figure 13 . False-color composite generated using topographic openness (green: 51×51 kernel), topographic openness (blue: 501×501 kernel) and slope (red: 3×3 kernel). Reddish colors represent high slope and low openness and greenish/bluish colors depict low slopes and relatively high openness (e.g., in flat and open areas).

Even though the field of spatial databases is more than 40 years old, most existing logical data models are highly focused either on spatial objects (vector data models) or spatial fields (raster data models). Furthermore, spatial index structures and query algorithms are still proposed for one of the approaches and little research work has been dedicated to index structures and query algorithms where both types of information are needed. However, due to the current high availability of different types of data, it is much more common nowadays that applications require querying vector and raster data at the same time.

This paper presents a method to perform a spatial query between a vector data set represented using an R-tree and a raster data set represented using a compact and space-efficient data structure called k 2 -tree that saves main memory space. Therefore, the method described in this paper solves two problems: first, it can be used to evaluate queries between vector and raster data without having to convert one of the data sets to the other data model and second, it saves main memory space, thus obtaining a more scalable system.


Raster* object. Typically a multi-layer type (RasterStack or RasterBrick)

fitted model of any class that has a 'predict' method (or for which you can supply a similar method as fun argument. E.g. glm, gam, or randomForest

character. Optional output filename

function. Default value is 'predict', but can be replaced with e.g. (depending on the type of model), or your own custom function.

Extent object to limit the prediction to a sub-region of x

data.frame. Can be used to add a constant for which there is no Raster object for model predictions. Particularly useful if the constant is a character-like factor value for which it is currently not possible to make a RasterLayer

integer. To select the column(s) to use if predict.'model' returns a matrix with multiple columns

logical. Remove cells with NA values in the predictors before solving the model (and return a NA value for those cells). This option prevents errors with models that cannot handle NA values. In most other cases this will not affect the output. An exception is when predicting with a boosted regression trees model because these return predicted values even if some (or all!) variables are NA

Lab Assignment

Please use the assignment page to complete your assignment and upload it to Stellar.

Created by Raj Singh and Joseph Ferreira.
Modified for 1999-2008 by Joseph Ferreira, Thomas H. Grayson, Jeeseong Chung , Jinhua Zhao, Xiongjiu Liao, Diao Mi, Michael Flaxman, and Yang Chen.
Last Modified: October 24, 2010 by Joe Ferreira.

Back to the 11.520 Home Page .
Back to the CRON Home Page.


While many operations in professional GIS’s are devoted to
producing professional looking maps, basic R is devoted
more to analytical operations.
A range of mathematical operations can be applied to matrices and vectors
to help to prepare spatial data and answer statistical questions.


The first operation is rasterizing — plotting points onto a raster.
This is exceedingly easy using the indexing operations in R
as shown below.

Histograms are an essential construct for examining the
distributions of values. In the histogram of classes in a raster layer.
a large proportion of values are 255. This is due to using the value
of 255 to represent ocean.


In ecological niche modeling of spatial data it
is usually necessary to mask out values that are irrelevant to the analysis.
One way to do it is to set a specific value in a mask vector, such as zero, and use
arithmetic operations to nullify the values in another vector, such as multiplication
by zero.

When using image processing operations, masking can be achieved
with addition or subtraction operations.

For example, in performing a masking operation, making use of the limited
range of the single byte in each cell with the following command, if the areas
to be masked have value 255 (white) then all areas in the image to be masked,
back.ppm will have the value 255 after the operation.


Proximity is often an important relationship to capture in ecological
analysis. Convolution is the operation whereby we use neighboring values to
determine the value of the central cell. In R a filter can be applied to
achieve convolution in one direction. The matrix can then
be transposed to apply the filter in the other directions.

Image processing packages usually have a smoothing operation
that achieves the same purpose. The program in neppbm is pnmsmooth
or pnmconvol. The utility pnmconvol uses a convolution matrix file.
In netpbm this is specified
with a pgm image as follows. The parameters width and height specify the
dimensions of hte convolution matrix.


Cropping refers to the trimming of unwanted parts of a 2D matrix to
leave the parts necessary for analysis. Foe example, a continental map may be
cropped to remove the area surrounding the set of points where a species occurs, so
as to develop more specific regional models.

This is done by pamcut in netpbm as follows:

The approach in R is to develop an array of indices for each of the cells in
the rectangle required with a command called expand.grid. This is however
not suitable for large data matrices.

Generalization or model development

All models are essentially generalizations, or simplification that enable expression
of theories in mathematical terms. As such, one of the simplest forms of generalization
is categorization, or clustering, where a large number of dissimilar items are sorted into a smaller
number of bins, based on their similarity. Once a set of bins or categories is established,
and there is a basis for deciding into which bin new items should go, new items can
be categorized. In this way, a categorization, or clustering can serve as a predictive

In R the basic operation for clustering is called kmeans. In kmeans, the data to be clustered
is partitioned into k groups such that the sum
of squares from points to the assigned cluster centres is minimized.
At the minimum, all cluster centres are at the mean of
the set of data points with the same category.

A similar operation is used in image processing called color quantization
or color reduction. Reducing the number of colors will compress the size of
an image at the expense of the number of colors. The utility in netpbm
for this purpose is called ppmquant.


Prediction can be achieved very efficiently with index mapping where
the dependent variable is one dimensional.
Index mapping is another fundamental operation for
mapping the values of a matrix into another set of values.
This operation can be used to apply the results of a model,
for example, changing the values in the cells from their
original value, to a probability for each specific class.

In image processing, index mapping is called palette mapping,
a process where the original colors are substituted for a new set.
Here a palette file defining the mapping is used.

The utility in the netpbm package is pamlookup, invoked with an image as
a lookup table for mapping the old colors to the new.

Note that this is a very efficient operation as the data in the image
does not change, only the small set of values in the palette of the image.

Joseph AE, Phillips DR: Accessibility and Utilization: Geographical Perspectives on Health Care Delivery. 1984, Harper & Row Ltd, London

Arcury TA, Gesler WM, Preisser JS, Sherman J, Spencer J, Perin J: The Effects of Geography and Spatial Behavior on Health Care Utilization among the Residents of a Rural Region. Health Services Research. 2005, 40: 135-156. 10.1111/j.1475-6773.2005.00346.x.

Cromley EK, McLafferty S: GIS and Public Health. 2002, Guilford Press, New York

Guagliardo M: Spatial accessibility of primary care: concepts, methods and challenges. International Journal of Health Geographics. 2004, 3 (3): 1-13.

McGuirk MA, Porell FW: Spatial patterns of hospital utilization: the impact of distance and time. Inquiry. 1984, 21: 84-95.

Shannon GW, Skinner JL, Bashshur RL: Time and distance: the journey for medical care. International Journal of Health Services. 1973, 3 (2): 237-244. 10.2190/FK1K-H8L9-J008-GW65.

McLafferty SL: GIS and Health Care. Annual Reviews in Public Health. 2003, 24: 25-42. 10.1146/annurev.publhealth.24.012902.141012.

Apparicio P, Abdelmajid M, Riva M, Shearmur R: Comparing alternative approaches to measuring the geographical accessibility of urban health services: Distance types and aggregation-error issues. International Journal of Health Geographics. 2008, 7 (7): 1-14.

Haynes R, Jones A, Sauerzapf V, Zhao H: Validation of travel times to hospital estimated by GIS. International Journal of Health Geographics. 2006, 5 (40): 1-8.

Phibbs CS, Luft HS: Correlation of Travel Time on Roads versus Straight Line Distance. Medical Care Research and Review. 1995, 52 (4): 532-542. 10.1177/107755879505200406.

Jones SG, Ashby AJ, Momin SR, Naidoo A: Spatial Implications Associated with Using Euclidean Distance Measurements and Geographic Centroid Imputation in Health Care Research. Health Services Research. 2010, 45: 316-327. 10.1111/j.1475-6773.2009.01044.x.

Couclelis H: People Manipulate Objects (but Cultivate Fields): Beyond the Raster-Vector Debate in GIS. Theories and Methods of Spatio-Temporal Reasoning in Geographic Space, Volume 639 of Lecture Notes in Computer Science. Edited by: Frank AV, Campari I, Formentini U. 1992, Springer, Berlin / Heidelberg, 65-77.

van Bemmelen J, Quak W, van Hekken M, van Oosterom P: Vector vs. Raster-based Algorithms for Cross Country Movement Planning. Proceedings Auto-Carto, Volume 11. 1993, 304-317.

Goodchild MF, Yuan M, Cova TJ: Towards a general theory of geographic representation in GIS. International Journal of Geographical Information Science. 2007, 21 (3): 239-260. 10.1080/13658810600965271.

Aday LA, Andersen R: A Framework tor the Study of Access to Medical Care. Health Services Research. 1974, 9 (3): 208-220.

Penchansky R, Thomas JW: The Concept of Access: Definition and Relationship to Consumer Satisfaction. Medical Care. 1981, 19 (2): 127-140. 10.1097/00005650-198102000-00001.

Khan AA: An integrated approach to measuring potential spatial access to health care services. Socio-Economic Planning Sciences. 1992, 26 (4): 275-287. 10.1016/0038-0121(92)90004-O.

Higgs G: A Literature Review of the Use of GIS-Based Measures of Access to Health Care Services. Health Services and Outcomes Research Methodology. 2004, 5 (2): 119-139. 10.1007/s10742-005-4304-7.

Higgs G: The role of GIS for health utilization studies: literature review. Health Services and Outcomes Research Methodology. 2009, 9 (2): 84-99. 10.1007/s10742-009-0046-2.

Martin D, Wrigley H, Barnett S, Roderick P: Increasing the sophistication of access measurement in a rural healthcare study. Health & Place. 2002, 8: 3-13. 10.1016/S1353-8292(01)00031-4.

Pedigo AS, Odoi A: Investigation of Disparities in Geographic Accessibility to Emergency Stroke and Myocardial Infarction Care in East Tennessee Using Geographic Information Systems and Network Analysis. Annals of Epidemiology. 2010, 20 (12): 924-930. 10.1016/j.annepidem.2010.06.013.

Shahid R, Bertazzon S, Knudtson M, Ghali W: Comparison of distance measures in spatial analytical modeling for health service planning. BMC Health Services Research. 2009, 9 (200): 1-14.

Witlox F: Evaluating the reliability of reported distance data in urban travel behaviour analysis. Journal of Transport Geography. 2007, 15 (3): 172-183. 10.1016/j.jtrangeo.2006.02.012.

Longley PA, Goodchild M, Maguire DJ, Rhind DW: Geographic Information Systems and Science. 2010, Wiley, Chichester

Dai D: Black residential segregation, disparities in spatial access to health care facilities, and late-stage breast cancer diagnosis in metropolitan Detroit. Health & Place. 2010, 16 (5): 1038-1052. 10.1016/j.healthplace.2010.06.012.

Schuurman N, Berube M, Crooks VA: Measuring potential spatial access to primary health care physicians using a modified gravity model. Canadian Geographer. 2010, 54: 29-45. 10.1111/j.1541-0064.2009.00301.x.

Wan N, Zhan FB, Zou B, Chow E: A relative spatial access assessment approach for analyzing potential spatial access to colorectal cancer services in Texas. Applied Geography. 2011, 32 (2): 291-299.

Kwan MP, Hong XD: Network-Based Constraints-Oriented Choice Set Formation Using GIS. Geographical Systems. 1998, 5: 139-162.

Lopez-Quilez A, Munoz F: Geostatistical computing of acoustic maps in the presence of barriers. Mathematical and Computer Modelling. 2009, 50 (5-6): 929-938. 10.1016/j.mcm.2009.05.021.

Messina JP, Shortridge AM, Groop RE, Varnakovida P, Finn MJ: Evaluating Michigan’s community hospital access: spatial methods for decision support. International Journal of Health Geographics. 2006, 5: 1-18. 10.1186/1476-072X-5-1.

Ray N, Ebener S: AccessMod 3.0: computing geographic coverage and accessibility to health care services using anisotropic movement of patients. International Journal of Health Geographics. 2008, 7: 1-17. 10.1186/1476-072X-7-1.

Tanser F, Gijsbertsen B, Herbst K: Modelling and understanding primary health care accessibility and utilization in rural South Africa: An exploration using a geographical information system. Social Science & Medicine. 2006, 63 (3): 691-705. 10.1016/j.socscimed.2006.01.015.

Sander HA, Ghosh D, van Riper D, Manson SM: How do you measure distance in spatial models? An example using open-space valuation. Environment and Planning B: Planning and Design. 2010, 37 (5): 874-894. 10.1068/b35126.

Upchurch C, Kuby M, Zoldak M, Barranda A: Using GIS to generate mutually exclusive service areas linking travel on and off a network. Journal of Transport Geography. 2004, 12: 23-33. 10.1016/j.jtrangeo.2003.10.001.

Michigan Department of Community Health: Certificate of Need Review Standards for Hospital Beds. 2009, [],

Birkmeyer JD, Siewers AE, Marth NJ, Goodman DC: Regionalization of High-Risk Surgery and Implications for Patient Travel Times. The Journal of the American Medical Association. 2003, 290 (20): 2703-2708. 10.1001/jama.290.20.2703.

Nallamothu BK, Bates ER, Wang Y, Bradley EH, Krumholz HM: Driving Times and Distances to Hospitals With Percutaneous Coronary Intervention in the United States: Implications for Prehospital Triage of Patients With ST-Elevation Myocardial Infarction. Circulation. 2006, 113 (9): 1189-1195. 10.1161/CIRCULATIONAHA.105.596346.

Berke E, Shi X: Computing travel time when the exact address is unknown: a comparison of point and polygon ZIP code approximation methods. International Journal of Health Geographics. 2009, 8 (23): 1-9.

Michigan Office of Highway Safety Planning: Establishing Realistic Speed Limits. Tech. Rep. OHSP 894, Lansing, MI

Wang F, Xu Y: Estimating O–D travel time matrix by Google Maps API: implementation, advantages, and implications. Annals of GIS. 2011, 17 (4): 199-209.

Price M: Slopes, Sharp Turns, and Speed. ArcUser. 2008,, (Spring):50–55

Price M: Convincing the Chief. ArcUser. 2009,, (Spring):50–54

Schuurman N, Fiedler R, Grzybowski S, Grund D: Defining rational hospital catchments for non-urban areas based on travel-time. International Journal of Health Geographics. 2006, 5 (43): 1-8.

Current JR, Schilling DA: Analysis of Errors Due to Demand Data Aggregation in the Set Covering and Maximal Covering Location Problems. Geographical Analysis. 1990, 22 (2): 116-126.

Watch the video: Reclassify, Overlay and Map Algebra of Raster Data in ArcGIS (October 2021).