Introduction

The AQFx model (Air Quality Forecast) is a chemical transport model that is run daily, over a range of nested domains covering Australia, or parts thereof. Here we use the SE, Syd and VTas domains, covering the whole continent, the Sydney region, and the Victorian/Tasmanian region. We start by expanding each grid to the largest extent, resampling all the grids to the same resolution using bilinear interpolation (0.3 degrees), and then overlaying them, so the finer scale Syd and VTas domains sit on top of the national domain.

aus <- read_sf("E:/geodata/Aus_Coastline/australia/aust.gpkg")
orast <- raster("data/current_merge.tif")
map <- tm_shape(orast) + tm_raster(breaks=c(-Inf,10,25,100,250,1000,Inf), palette=c("#4CAF50FF","#FFEB3BFF","#FF9800FF","#F44336FF","#673AB7FF","#000000FF"),title="AQFx") + tm_shape(aus) + tm_borders()
print(map)

We can see significant smoke concentration across northern Australia, including severe levels in south-east Queensland.

The next step is to read in our station data, and calcuate differences between monitor PM and model PM for the current hour.

dat <- read_csv("data/dat.csv")

# Extract model PM
dat$mod_pm <- extract(orast,cbind(dat$Lon,dat$Lat))

# Calculate differences (residuals)
dat$resid <- dat$mod_pm - dat$PM25

# Convert to sp (spatial) object, needed for gstat modeling
coordinates(dat) <- ~Lon+Lat
projection(dat) <- projection(orast)

We will try two different spatial models, one using a krieg, and one using simple IDW.

# Krieg of the residuals - untransformed
m <- autofitVariogram(resid~1,dat)
mod_k <- gstat(NULL, "p_resid", resid~1, dat, model=m$var_model)

# IDW  the residuals - untransformed
mod_idw <- gstat(NULL, "p_resid", resid~1, dat,nmax=15)


# Interpolate them both across our raster

i_krieg <- interpolate(orast,mod_k)

## [using ordinary kriging]

i_idw <- interpolate(orast,mod_idw)

## [inverse distance weighted interpolation]

# Adjust raster by subtracting the smooth

orast_adj <- orast - i_krieg
orast_adj_2 <- orast - i_idw

# We can get values below zero - clip them

values(orast_adj)[values(orast_adj)<0]=0
values(orast_adj_2)[values(orast_adj_2)<0]=0

If we compare the krieg and IDW interpolations, we can see the krieg is smoother, while IDW is tighter to local variation. In theory, the krieg should give the best interpolation, as it takes into account spatial autocorrelation and structure. However, we are using an autokrieg algorithm to calcualte the variogram, rather than doing it in an informed way with a specific spatial scale and structure.

m1 <- tm_shape(i_krieg) + tm_raster(title = "Krieg") + tm_shape(aus) + tm_borders()
m2 <- tm_shape(i_idw) + tm_raster(title = "IDW") + tm_shape(aus) + tm_borders()
mc <- tmap_arrange(m1,m2,ncol=2)
print(mc)

Now we plot the input AQFx surface by these adjustments; we can see the big smoke plume across northern Australia is “pulled down” towards reality, although not completely - it remains in areas where there are few monitoring stations.

m1 <- tm_shape(orast_adj) + tm_raster(breaks=c(-Inf,10,25,100,250,1000,Inf), palette=c("#4CAF50FF","#FFEB3BFF","#FF9800FF","#F44336FF","#673AB7FF","#000000FF"),title="Krieg") + tm_shape(aus) + tm_borders()
m2 <- tm_shape(orast_adj_2) + tm_raster(breaks=c(-Inf,10,25,100,250,1000,Inf), palette=c("#4CAF50FF","#FFEB3BFF","#FF9800FF","#F44336FF","#673AB7FF","#000000FF"),title="IDW") + tm_shape(aus) + tm_borders()
mc <- tmap_arrange(m1,m2,ncol=2)
print(mc)

We want to ensure the areas exactly around monitoring stations (within 20m) exactly reflect the values of the monitoring stations. To do this, we create a mask with 20km areas around each monitor, interpolate across the raw monitor values, and then assign the areas within those masks with the raw monitor interpolation.

# Create 20km mask around stations
br <- buffer(dat,20000)
tmp <- orast
values(tmp)<-1
tmp <- mask(tmp,br)
values(tmp)[is.na(values(tmp))]<-0
dfp2 <- tmp
dfp2<-1-tmp

# Take the krieg-adjusted output, mask out holes
orast_adj <- orast_adj * dfp2


# Do IDW-interpolation of PM2.5
# Because we are doing raw values, not residuals, log the data
dat$lPM <- log1p(dat$PM25)

gOK2 <- gstat(NULL, "lPM", lPM~1, dat, nmax=15, set=list(idp = 1.5))
OK2 <- raster::interpolate(orast_adj, gOK2)

## [inverse distance weighted interpolation]

# Exponentiate to return to normal values
OK2<-expm1(OK2)

# Fill the holes in the mask with these interpolated station values
OK<-OK2*tmp

# Add the two rasters - the adjustedmodel and the masked interpolation - together
OK <- OK + orast_adj

Finally we can plot the adjusted model with the station values placed around stations.

map <- tm_shape(OK) + tm_raster(breaks=c(-Inf,10,25,100,250,1000,Inf), palette=c("#4CAF50FF","#FFEB3BFF","#FF9800FF","#F44336FF","#673AB7FF","#000000FF"),title="Final Output") + tm_shape(aus) + tm_borders()
print(map)