[R-sig-Geo] Raster processing, pop density estimation and more

[Luciano La Sala]

Dear r-sig-geo people,


I have a raster image of human population counts; a LandScan grid at a
resolution of 1 km where each pixel's value represents population number,
which I am processing using Quantum GIS and GRASS the modules inside of it.


Using R, I need to: 

1) Take each cell of the raster one by one and retrieve the population count
for that cell.

2) Work out the number of households per cell: for this, I should take a
random draw from a Poisson distribution assuming that there is a mean of 4
people per household. Keep a running total individuals allocated to
households so that when we reach the total population size for the given
cell the process stops. The total number of Poisson draws taken for a cell
will equal the number of households.

3) Finally, I am interested in the number of households in a cell that keep
poultry. I should take a random draw from a binomial distribution. For this
the number of trials would be the total number of households and the
probability of keeping poultry would depend on human population density. If
population density is high (i.e. the household is in a city) the probability
of keeping poultry is zero. 

I colleague who did a very similar work (only that he worked with ArcView)
wrote the following code below to get the job done:  

----------------------------------------------------------------------------
---------------------------------------------------

# Estimating the number of poultry-positive parcels per grid cell

# Create a grid theme of population density - based on the grid themes
specifying population counts and area.
# Export population count and population density files from ArcView to *.asc
files.

pop.file <- "vc_pop.asc"
den.file <- "vc_den.asc"
out.file <- "vc_byp.asc"

out.xllcorner <- -74.010559082031
out.yllcorner <- -33.743896484375
out.cellsize <- 0.008333333333333
out.na <- -9999

dat.pop <- as.matrix(read.table(file = pop.file, sep = " ", na.strings =
"-9999", skip = 6)) dat.den <- as.matrix(read.table(file = den.file, sep = "
", na.strings = "-9999", skip = 6))

# Set up a matrix to collect results:
dat.byp <- dat.pop
dat.byp[,][dat.byp[,] >= 0] <- 0

# 2. Loop through each grid cell of population file and work out the
estimated number of establishments in each:
nrows <- dim(dat.den)[1]
ncols <- dim(dat.den)[2]

for(i in 1:nrows){
  for(j in 1:ncols){
    total <- 0

    # Only proceed if data is present:
    if(!is.na(dat.pop[i,j]) & !is.na(dat.den[i,j]) ){

      # What is the population size and density for this cell?
      tmp.pop <- dat.pop[i,j]
      tmp.den <- dat.den[i,j]

      # What is the population density category for this cell?
      id.den <- 0
      if(tmp.den < 1.50) id.den <- 1
      if(tmp.den >= 1.50 & tmp.den < 4.20) id.den <- 2
      if(tmp.den >= 4.20 & tmp.den < 6.00) id.den <- 3
      if(tmp.den >= 6.00 & tmp.den < 10.00) id.den <- 4
      if(tmp.den >= 10.00) id.den <- 5

      while(sum(total) < tmp.pop) {
       n <- rpois(n = 1, lambda = 4)
       # Vector of the number of individuals in each establishment:
       total <- c(total, n)
       }

      # Number of establishments in this cell:
      n.estab <- length(total) - 1

      # Estimated number of establishments that keep poultry (set to vary
according to human population density):
      if(id.den == 1){dat.byp[i,j] <- rbinom(n = 1, size = n.estab, prob =
0.90)}
      if(id.den == 2){dat.byp[i,j] <- rbinom(n = 1, size = n.estab, prob =
0.50)}
      if(id.den == 3){dat.byp[i,j] <- rbinom(n = 1, size = n.estab, prob =
0.20)}
      if(id.den == 4){dat.byp[i,j] <- rbinom(n = 1, size = n.estab, prob =
0.10)}
      if(id.den == 5){dat.byp[i,j] <- 0}

      # if(id.den == 1){dat.byp[i,j] <- qbinom(0.50, size = n.estab, prob =
0.90)}
      # if(id.den == 2){dat.byp[i,j] <- qbinom(0.50, size = n.estab, prob =
0.50)}
      # if(id.den == 3){dat.byp[i,j] <- qbinom(0.50, size = n.estab, prob =
0.05)}
      # if(id.den == 4){dat.byp[i,j] <- 0}
    }
   }
  }

rval.byp <- dat.byp

----------------------------------------------------------------------------
---------------------------------------------------

I am using Quantum GIS/GRASS and since my original raster (Landscan grid) is
at resolution 1 km, then the population count and density (per sq km) should
be the same. Using the GRASS module r.out.xyz I got an ASCII text file with
three columns: X, Y, value. The X-Y are cell enters, and value would
represent Population count.
 
That being said, I guess that my pop density grid would be the LandScan grid
as it is (people/sq km). If this is correct, then I see no point in
following the steps below (as described above)?   

# Create a grid theme of population density - based on the grid themes
specifying population counts and area.
# Export population count and population density files from ArcView to *.asc
files.

Any help with this will be greatly appreciated! 

Luciano La Sala