¿Entender la resolución al rasterizar en R?

Question

Me estoy metiendo en la construcción de superficies - útil para una serie de proyectos de salud pública - y parece que no puedo dar sentido a lo que está sucediendo con la resolución y el proceso de construcción de raster para conseguir lo que quiero. Tengo un ejemplo que espero que sea trivial y reproducible a continuación para construir. Por supuesto, por razones de proyección, esto podría no ser sensible a esta escala (lo que realmente estoy haciendo es más en el bloque / grupo de bloques / nivel de contrato), pero creo que esto hace que para un punto de discusión digerible.

Me gustaría tener la población dentro de un gran círculo en los Estados Unidos. Para ello, tengo un shapefile, algunos datos de población y un círculo de unos 1000 km en el centro del país.

¿Cómo puedo calcular la población en ese círculo?

Imagen y código roto abajo.

¿Qué opinas?

  #################################################################
  # "I don't get rasters"
  # an example by Mike Dolan Fliss. NC Epidemiology.
  #################################################################

  #NOTE: You'll need to install.packages() these if you don't use them.
  library(rgdal)
  library(sp)
  library(dplyr)
  library(ggplot2) # Could get better plots with ggplot, but for speed we'll use base/sp
  library(rgeos)
  library(raster)
  library(xlsx) #For the description table, an excel file

  # Download, read, and touch up the geographic data ##############
  #################################################################
  download.file("http://www2.census.gov/geo/tiger/TIGER2010DP1/State_2010Census_DP1.zip", "states.zip") #~9meg
  unzip("states.zip")
  states.spdf = readOGR(dsn = ".", layer = "State_2010Census_DP1", stringsAsFactors = F)
  states.spdf = states.spdf[!(states.spdf$NAME10 %in% c("Alaska", "Hawaii", "Puerto Rico")), ] #get continental

  descriptions.df = read.xlsx("DP_TableDescriptions.xls", 1)
  head(descriptions.df) # Total pop is in : DP0010001
  states.spdf$total.pop = states.spdf$DP0010001 #For readability

  plot(states.spdf) #Plot in base
  spplot(states.spdf, "total.pop") #Plot in sp, gorgeous Lisa Frank colors
  spplot(states.spdf, "total.pop", col.regions=rev(heat.colors(20))) #heat

  # Project & build some buffers ##################################
  #################################################################
  epsg_codes = make_EPSG()# Let's go Albers equal area. I'm not great with projections.  
  us.equalarea.prj = epsg_codes$prj4[epsg_codes$note=="# US National Atlas Equal Area"]
  states.spdf = spTransform(states.spdf, CRS(us.equalarea.prj)) #unit is now m.

  center = gCentroid(states.spdf)
  little.buffer.spdf = gBuffer(center, width=1000*1000) #1000km radius "circle"
  big.buffer.spdf = gBuffer(center, width=2700*1000) #2700 radius - all of US
  plot(states.spdf);plot(big.buffer.spdf, add=T); plot(little.buffer.spdf, add=T)

  # Now, I fail to make the right raster... :) ####################
  #################################################################
  pop.raster = raster(extent(states.spdf))
  e = extent(states.spdf); (e@xmax-e@xmin)/1000 #US extent is 4560km across. ok...
  projection(pop.raster) = proj4string(states.spdf)
  raster.res = 1000*100 #makes 100km x 100km raster squares (?)
  res(pop.raster) = c(raster.res, raster.res) 
  states.spdf$frac.pop = states.spdf$total.pop / (states.spdf$ALAND10+states.spdf$AWATER10)*raster.res
  #head(bgs.spdf$frac.pop)
  pop.raster = rasterize(states.spdf, pop.raster, "frac.pop")

  plot(pop.raster) #Plot - looks good
  plot(states.spdf, border="black", add=T)

  cellStats(pop.raster, sum) #3093? Clearly I'm confused
  sum(states.spdf$total.pop) #US pop = 307 million

  # ... and subsequently fail to get the right area ###############
  #################################################################
  plot(rp)
  plot(states.spdf, border="black", add=T)
  plot(little.buffer.spdf, border="blue", add=T)
  # What's the area in this blue circle?

  e = extract(pop.raster, big.buffer.spdf)
  sum(unlist(e), na.rm=T) #3093... somethings
  #^ Here we are, back to the same as cellStats sum, above.

  e = extract(pop.raster, little.buffer.spdf)
  sum(unlist(e)) #663... 
  #^ I'd like this to represent the population contained in the circle

Answer 1

Hay dos problemas con su código.

La primera y más crucial se muestra a continuación:

states.spdf$frac.pop = states.spdf$total.pop / (states.spdf$ALAND10+states.spdf$AWATER10)*raster.res

Estás asumiendo que la población se distribuye por igual en el espacio, sin embargo tratas la densidad de forma errónea. En lugar de utilizar raster.res para distribuir la población raster.res^2 ya que se trata de un área. La línea de código correcta se da a continuación:

states.spdf$frac.pop = states.spdf$total.pop / (states.spdf$ALAND10+states.spdf$AWATER10)*raster.res^2

Implementar esta línea te acercaría al número real:

sum(states.spdf$total.pop) #Porcentaje de población estadounidense = 307 millones

[1] 306675006

cellStats(pop.raster, sum)

[1] 309351269

El segundo problema está relacionado con la resolución en su efecto sobre los límites:

Una célula de 100 KM probablemente superará dos estados. La reducción del tamaño de la celda también disminuye la aparición de estas superposiciones. Por lo tanto, mejora la asignación de población a las celdas a lo largo de los límites. Véanse los resultados para 10 KM a continuación:

cellStats(pop.raster, sum) # 10 KM de tamaño de celda

[1] 305963050

sum(states.spdf$total.pop) #Porcentaje de población estadounidense = 307 millones

[1] 306675006

Además, tiene un efecto significativo en los resultados visuales.

Por último, hay que tener en cuenta que existe un compromiso entre la resolución y el tiempo de ejecución, por lo que la resolución debe elegirse con cuidado.