I am having 2 sets of raster
data and their names are:
-
ntl_'a number'.tif
-
pop_'a number'.tif
My goal is to create a function that reads the filenames (e.g., ntl_1.tif and pop_1.tif) and then perform the below analysis to raster with the same number in their layer names:
library(raster)
library(DescTools)
#create a data.frame of values from the NTL and pop raster data
ntl = raster("path/ntl_1.tif")
vals_ntl <- as.data.frame(values(ntl))
ntl_coords = as.data.frame(xyFromCell(ntl, 1:ncell(ntl)))
combine <- as.data.frame(cbind(ntl_coords,vals_ntl))
pop<-raster("path/pop_1.tif")
pop = resample(pop, ntl, method = 'bilinear')
vals_pop <- as.data.frame(values(pop))
block.data <- as.data.frame(cbind(combine, vals_pop))
names(block.data)[3] <- "ntl"
names(block.data)[4] <- "pop"
block.data <- na.omit(block.data)
block.data = subset(block.data, select = -c(x, y))
# sort by ntl
block.data <-block.data[order(block.data$ntl),]
ntl_vector <- block.data[ , "ntl"]
pop_vector <- block.data[ , "pop"]
#compute gini index
Gini(ntl_vector, pop_vector, unbiased = FALSE)
The above code is when using a pair of raster. In my case I have hundreds of pairs so I'd like to automate the process and hopefully to get the results (i.e., the gini coefficient) of every pair in my console or, even better, in a data.frame
.
I made a lot of tries to create that function so I don't know which of my tries should I post. Here are some dummy data:
new("RasterStack", filename = "", layers = list(new("RasterLayer",
file = new(".RasterFile", name = "", datanotation = "FLT4S",
byteorder = "little", nodatavalue = -Inf, NAchanged = FALSE,
nbands = 1L, bandorder = "BIL", offset = 0L, toptobottom = TRUE,
blockrows = 0L, blockcols = 0L, driver = "", open = FALSE),
data = new(".SingleLayerData", values = 1:9, offset = 0,
gain = 1, inmemory = TRUE, fromdisk = FALSE, isfactor = FALSE,
attributes = list(), haveminmax = TRUE, min = 1L, max = 9L,
band = 1L, unit = "", names = "layer.1"), legend = new(".RasterLegend",
type = character(0), values = logical(0), color = logical(0),
names = logical(0), colortable = logical(0)), title = character(0),
extent = new("Extent", xmin = -180, xmax = 180, ymin = -90,
ymax = 90), rotated = FALSE, rotation = new(".Rotation",
geotrans = numeric(0), transfun = function ()
NULL), ncols = 3L, nrows = 3L, crs = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs"),
history = list(), z = list()), new("RasterLayer", file = new(".RasterFile",
name = "", datanotation = "FLT4S", byteorder = "little",
nodatavalue = -Inf, NAchanged = FALSE, nbands = 1L, bandorder = "BIL",
offset = 0L, toptobottom = TRUE, blockrows = 0L, blockcols = 0L,
driver = "", open = FALSE), data = new(".SingleLayerData",
values = 1:9, offset = 0, gain = 1, inmemory = TRUE, fromdisk = FALSE,
isfactor = FALSE, attributes = list(), haveminmax = TRUE,
min = 1L, max = 9L, band = 1L, unit = "", names = "layer.2"),
legend = new(".RasterLegend", type = character(0), values = logical(0),
color = logical(0), names = logical(0), colortable = logical(0)),
title = character(0), extent = new("Extent", xmin = -180,
xmax = 180, ymin = -90, ymax = 90), rotated = FALSE,
rotation = new(".Rotation", geotrans = numeric(0), transfun = function ()
NULL), ncols = 3L, nrows = 3L, crs = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs"),
history = list(), z = list())), title = character(0), extent = new("Extent",
xmin = -180, xmax = 180, ymin = -90, ymax = 90), rotated = FALSE,
rotation = new(".Rotation", geotrans = numeric(0), transfun = function ()
NULL), ncols = 3L, nrows = 3L, crs = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs"),
history = list(), z = list())