0

I want to calculate resistance distance based on circuit theory for a set of species for the city of Berlin, GER. For the city area I use land cover, sealed soil and building height data. This data is freely available via http://fbinter.stadt-berlin.de/fb/index.jsp. I use ArcMap 10.3.
As proposed by Koen et al. 2010 (http://dx.doi.org/10.1371/journal.pone.0011785 ). To sum up the paper, I want to create a buffer around Berlin to account for edge effects. This buffer should show the same frequency distribution of resistance values as the city raster.

For example:
Scenario - Frequency Distribution (low/medium/high)

Berlin - (46/50/98)
Random - (46/50/98)

I started an attempt using R 3.3.3 on MacOS 10.11.6. So far I can import the sealed.tif (which is the city area), transform to a matrix and get the frequency distribution. I then created a new randomized matrix with values from 0 to 100.

setwd("/Users/franzxaverpoellinger/Desktop/Statistics/Raster/raster")

library(raster) 
library(sp)
library(rgdal)
library(raster)

str_name<-'sealed.tif' 
sealed=raster(str_name)
plot(sealed)

# convert from .tif to a matrix
sealed.matrix <- as.matrix(sealed)

seal.cut = cut(sealed.matrix, breaks = seq(0, 100, by = 10), right=FALSE)
seal.cut.freq = table(seal.cut)
(seal.cut.freq/sum(seal.cut.freq))*100 # this is the aimed frequency

# and now the same for the buffer .tif
str_name.buffer <-'buffer.tif' 
buffer=raster(str_name.buffer)
buffer.matrix <- as.matrix(buffer)
str(buffer.matrix)

# create a matrix the same size as buffer.hist
new.buffer.matrix <- matrix(sample.int(100, size= 457*537, replace = TRUE), nrow = 457, ncol = 537)
new.buffer.matrix

buffer.cut = cut(new.buffer.matrix, breaks = seq(0, 100, by = 10), right=FALSE)
buffer.cut.freq = table(buffer.cut)
(buffer.cut.freq/sum(buffer.cut.freq))*100

Now I want the new.buffer.matrix to show the same frequency distribution of as does sealing.matrix. If I had the same distribution I would (somehow) go back from matrix to .tif and then just clip it to the extend that I need.

1 Answer 1

1

I have a definite answer now, from input .tif down to exporting a .tif file for further use in QGIS or ArcGIS.

setwd("")

library(raster)
library(rgdal)

input_tif <-'any_tif_should_run.tif' # loading .tif file
input_raster <- raster(input_tif)
buffer_area <- readOGR('buffer_area.shp') # load shapefile of buffer if randomized area should be clipped to a certain extent
# if not remove masking in line 68

#### convert from .tif to a matrix
{
 input_matrix <- as.matrix(input_raster)
 run.projection  <- projection(input_raster)
 rows <- nrow(input_matrix) # number of rows
 cols <- ncol(input_matrix) # number of columns
 length.input_matrix <- length(input_matrix) # number of values in the input matrix
 input_matrix.cut <- cut(input_matrix, breaks = seq(0, 100, by = 10), right=FALSE) # cuts the matrix values by 10
 input_matrix.cut.freq <- table(input_matrix.cut) 
 str.input_matrix.cut.freq <- (input_matrix.cut.freq/sum(input_matrix.cut.freq))*100 # calculate relative frequency of values

# this produces the vector to fill the output.raster
 vector_10 <- c(1,11,21,31,41,51,61,71,81,91) # create vector
 input_matrix.freq.absolut <- round((str.input_matrix.cut.freq/100)*(length.input_matrix))
 input_vector <- rep(vector_10, input_matrix.freq.absolut) # vector with the values of
 # vector.10 with the absolute frequency of values of str.seal.cut.freq.absolut
 input_vector.sampled <- sample(input_vector) # randomize the vector values
}

# here the length of the vector to fill the raster with randomized values is cropped or extended
# to the number of values e.g. length in the matrix
if (length(input_vector.sampled) > length.input_matrix) {
 input_vector.sampled_out <- input_vector.sampled[1:length.input_matrix]
 } else if (length(input_vector.sampled) < length.input_matrix) {
 to.fill <- length.input_matrix-length(input_vector.sampled)
 appendix.vector <- rep(1, to.fill)
 input_vector.sampled.out <- c(input_vector.sampled, appendix.vector)
} else {
 input_vector.sampled.out <- input_vector.sampled
}
# create a matrix from input_vector.sampled.out in order to compare the frequency distribution with the input .tif/raster/matrix
{
approx.matrix <- matrix(input_vector.sampled.out,nrow = rows, ncol = cols) # create a matrix from vector to test relative frequency of classes
approx.cut <- cut(approx.matrix, breaks = seq(0, 100, by = 10), right=FALSE)
approx.cut.freq <- table(approx.cut)
str.input_vector <- (approx.cut.freq/sum(approx.cut.freq))*100
}
str.input_matrix.cut.freq    # compare distribution of input matrix and vector for raster conversion
str.input_vector

# convert matrix to raster
{
 x <- raster()
 ncol(x) <- cols
 nrow(x) <- rows
 projection(x) <- run.projection
 x <- setValues(x, input_vector.sampled.out) # fill with values of input_vector.sampled.out with random values
 new.x <- projectRaster(x, input_raster) # apply extent etc. of input raster
 new.x <- setValues(new.x, input_vector.sampled.out)
 covered.raster <- cover(input_raster, new.x) # cover NA's of 'sealed' with the random values of raster 'new.x'
 clipped <- mask(covered.raster, buffer_area) # clip covered.raster by buffer_area.shp
 plot(clipped)
}
writeRaster(clipped, filename = "output_raster.tif")

I wrote this script in a manner that it should run with any kind of .tif input.

3
  • Hey poelinf, Would you mind taking a look here? stackoverflow.com/questions/42680596/… Mar 17, 2017 at 15:45
  • Hi @VijayRamesh, I'm flattered by your request but have honestly no idea what you even try to do...never worked with that kind of stuff. Good luck!
    – blabbath
    Mar 19, 2017 at 11:57
  • No problem @poellinf :). Just thought I would ask. Looks like I figured it out! Mar 19, 2017 at 17:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.