In this tutorial we’ll learn the pixel level linear regression. Here, we have used black carbon (BC) and cloud fraction (CF) to run the regression. The BC is considered as independent variable while the CF is considered as dependent variable. The core idea is to learn the effect of BC on the cloud fraction. As BC is a radiation absorptive aerosols, a higher presence of BC in atmosphere leads to decrease in cloud fraction as they absorb the radiation leading to evaporate the cloud droplets. We’ll go by general expression of linear model as given below.

\(\hat{y_{i}} = \beta_{0} + \beta_{1}*x\)

The prime focus of the linear model is to map the \(\beta_{1}\) or the slope value which suggests the rate of change of \(y_{i}\) as \(x_{i}\) change by an unit. To execute the linear modelling we’ll use the time series BC and CF data from Giovanni webpage (https://giovanni.gsfc.nasa.gov/giovanni/). The CF was taken from MODIS, whereas BC was taken from MERRA-2 reanalysis data.

# import the seasonal BC data. Within '', the file address is given. this address may be different for your system. Make sure you double check it beforoe proceeding.
library(raster)
Loading required package: sp
library(rasterVis)
Loading required package: lattice
library(latticeExtra)

Attaching package: 'latticeExtra'
The following object is masked from 'package:rasterVis':

    horizonplot
mam_bc = brick('E:/humboldt_PD_data_work/bc_mass/mam_bc.nc')
Loading required namespace: ncdf4
jjas_bc = brick('E:/humboldt_PD_data_work/bc_mass/jjas_bc.nc')
on_bc = brick('E:/humboldt_PD_data_work/bc_mass/on_bc.nc')
djf_bc = brick('E:/humboldt_PD_data_work/bc_mass/djf_bc.nc')

# import the seasonal cloud fraction data
mam_cf = brick('E:/humboldt_PD_data_work/cloud_fract/mam_cf.nc')
jjas_cf = brick('E:/humboldt_PD_data_work/cloud_fract/jjas_cf.nc')
on_cf = brick('E:/humboldt_PD_data_work/cloud_fract/on_cf.nc')
djf_cf = brick('E:/humboldt_PD_data_work/cloud_fract/djf_cf.nc')

The data are now imported. Let’s stack them. The following few lines of code will stack all the data. Note that naming convention of the new file is in accordance with how they are been stacked. However, it is free to use as per the your convenience. Since this data is used over Indian region, the BC file should matched with the extension of CF file. Thus, crop function is used to make the extent of the two files same.

# stack the two seasonal time stacked images
# if different extent, use mask function
mam_bc_crop = crop(mam_bc,extent(mam_cf))

# now stack them
mam_cf_bc = stack(mam_cf,mam_bc_crop)

jjas_bc_crop = crop(jjas_bc,extent(jjas_cf))
jjas_cf_bc = stack(jjas_cf,jjas_bc_crop)

on_bc_crop = crop(on_bc,extent(on_cf))
on_cf_bc = stack(on_cf,on_bc_crop)

djf_bc_crop = crop(on_bc,extent(on_cf))
djf_cf_bc = stack(djf_cf,djf_bc_crop)

You can quickly check the details of every raster brick imported so far. Here is an example.

print(mam_cf_bc)
class      : RasterStack 
dimensions : 38, 37, 1406, 40  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 62, 99, 4, 42  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : X2001.04.01.1, X2002.04.01.1, X2003.04.01.1, X2004.04.01.1, X2005.04.01.1, X2006.04.01.1, X2007.04.01.1, X2008.04.01.1, X2009.04.01.1, X2010.04.01.1, X2011.04.01.1, X2012.04.01.1, X2013.04.01.1, X2014.04.01.1, X2015.04.01.1, ... 
min values :             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?, ... 
max values :             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?,             ?, ...

Let’s execute the regression function on the stacked image. We’ll write two functions; one for \(\beta_{1}\) and other for finding out the significance level (p-value) of the estimated \(\beta_{1}\).

# define function for slope
reg_fun = function(x){
  if(is.na(x[1])){NA}
  else{lm(x[1:20] ~ x[21:40])$coefficients[2]}
}

#fnction for p-value
reg_fun2 = function(x){
  if(is.na(x[1])){NA}
  else{m=lm(x[1:20] ~ x[21:40]);summary(m)$coefficients[8]}
}

# apply the function over the stacked image
sl_mam = calc(mam_cf_bc, reg_fun)
sl_jjas = calc(jjas_cf_bc, reg_fun)
sl_on = calc(on_cf_bc, reg_fun)
sl_djf = calc(djf_cf_bc, reg_fun)

# get p-value
p_mam = calc(mam_cf_bc, reg_fun2)
p_jjas = calc(jjas_cf_bc, reg_fun2)
p_on = calc(on_cf_bc, reg_fun2)
p_djf = calc(djf_cf_bc, reg_fun2)

Since we’re only interested to the p-value which is lower than 0.05 (p< 0.05), we need to remove all other p-values from the analysis. The purpose is to put significant marker on the estimated \(\beta_{1}\) which are statistically significant. Following few lines of code will first subset the p-value with a condition (p<0.05), then convert them into vector layer so that they could be put on the \(\beta_{1}\).

# define function to get p-value less than 0.05
f_sig = function(x){
  ifelse(x <= 0.05, x*1, NA)
}
p_mam_05 = calc(p_mam, f_sig)
p_jjas_05 = calc(p_jjas, f_sig)
p_on_05 = calc(p_on, f_sig)
p_djf_05 = calc(p_djf, f_sig)

# convert the p-val lt 0.05 into vector point
rtf_mam_05 = rasterToPoints(p_mam_05, spatial = T)
rtf_jjas_05 = rasterToPoints(p_jjas_05, spatial = T)
rtf_on_05 = rasterToPoints(p_on_05, spatial = T)
rtf_djf_05 = rasterToPoints(p_djf_05, spatial = T)

We’re almost close to our final outcome. Now, we need to stack the estimated slope together so that we can plot the slopes of all the seasons together.

slope_stack = stack(sl_mam,sl_jjas,sl_on,sl_djf)
# assign names
names(slope_stack) <- c('MAM','JJAS','ON','DJF')
library(sp)

# import the global admin
levelplot(slope_stack/1e6, layers = 1:4,
          colorkey = list(title = list('Slope',cex=0.8), space = 'right',
                                       height = 1, width = 1.3),
          par.settings = BuRdTheme())+layer(sp.points(rtf_mam_05, pch = 20, cex=0.5, col='black'), columns=1, rows = 1)+layer(sp.points(rtf_jjas_05,pch = 20, cex=0.5, col='black'), columns = 2, rows = 1)+layer(sp.points(rtf_on_05,pch = 20, cex=0.5, col='black'), columns = 1, rows = 2)+layer(sp.points(rtf_djf_05,pch = 20, cex=0.5, col='black'), columns = 2,rows = 2)