I am trying to get a raster file of nightlights data for Brazil and subdivide/regroup it into several rasters for each state, in a list or data.frame in R

I managed to download the nightlights file, import into R using the raster package, and crop it for a bounding box of Brazil. Reproducible code:


#Source: 'https://ngdc.noaa.gov/eog/data/web_data/v4composites/F182013.v4.tar'  
#Download this file manually ( untar() command does not work if I use download.file() to download)
fo <- 'C:/.../NightLights/'
f <- list.files(fo,full.names = T)
untar(f,exdir = fo)
tif <- paste0(fo,'/F182013.v4c_web.stable_lights.avg_vis.tif')
gunzip(paste0(tif,'.gz'),tif )

bra = getData('GADM', country="BRA", level=0) #country division
bra1 = getData('GADM', country="BRA", level=1) #state disions

F18_bra <- raster(tif) %>% crop(bra) #imports raster for a Bunding Box of Brazil

As the above code indicates I also imported the polygons for Brazilian states (object bra1)

I want to do a spatial join of the raster and the state polygons.

In the end, I want a data.frame, or list, with: state_id, state_raster. I.e., one line/element per state, with a raster object with the pixels inside the polygon of each state on each cell of state_raster (like a bucket or geom in Postgis)

Can I do this in R, preferably using the sp, sf and/or raster packages?

  • I suspect you want to rasterize() the polygon with some ID variable and that gives you a raster like F18_bra but with polygon ID values. You can then apply any function grouped over those ID values on the F18_bra raster. If you can make a reproducible example (maybe generate a fake "lights" raster of similar scale) and that doesn't use library(tidyverse) I'll write a full answer for you. – Spacedman May 11 at 7:47

you can use the combination of the crop() and mask() functions to crop the part of your raster that lies within each of the polygons of your bra1 object. Then, you might use lapply() to "loop over" all states to obtain a list, as you desired. Make sure to have both objects at the same CRS:

bra1 <- spTransform(bra1, tif@crs)
output <- lapply(1:nrow(bra1), function(x) {
  list(state_id = x, state_raster = mask(crop(tif,bra1[x,]), bra1[x,]))

The output, with a polygon and raster of mine, is a list as you desired:

> output
[1] 1

class       : RasterLayer 
dimensions  : 52, 75, 3900  (nrow, ncol, ncell)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -17.58333, -11.33333, 12.33333, 16.66667  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : act2000_r_mze_2000_har 
values      : 0, 0.3313416  (min, max)

[1] 2

class       : RasterLayer 
dimensions  : 151, 147, 22197  (nrow, ncol, ncell)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -17.08333, -4.833333, 14.75, 27.33333  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : act2000_r_mze_2000_har 
values      : 0, 0.2182628  (min, max)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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