knitr::opts_chunk$set(echo = TRUE) library(MazamaSatelliteUtils) library(MazamaSpatialUtils) setSatelliteDataDir("~/Data/Satellite") setSpatialDataDir("~/Data/Spatial") # get bbox for Texas loadSpatialData("USCensusStates") tx <- subset(USCensusStates, stateCode == "TX") bbox_tx <- sp::bbox(tx) goesaodc_downloadAOD(satID = "G16", datetime = "20191291201274", isJulian = TRUE, timezone = "UTC") nc <- goesaodc_openFile("OR_ABI-L2-AODC-M6_G16_s20191291201274_e20191291204047_c20191291210009.nc") pal_aod <- colorRampPalette(c("lightgoldenrod1", "red3")) pal_count <- colorRampPalette(c("powderblue", "dodgerblue4")) pal_sd <- colorRampPalette(c("lightpink", "darkred"))
A raster image divides geographic space into uniform pixels or "cells".
The uniform size of these cells is determined by the resolution of the raster so
that higher resolutions reduce cell area and increase the detail and accuracy of
the final image. The goesaodc_createRaster()
function defines a colored
cell for each grid space with at least one data point, which raises two issues:
When a grid space is sparsely populated, a cell will be defined for the whole area using only a few data possibly concentrated in just a small corner. These might not accurately represent the AOD level for the rest of the cell.
When a grid space is densely populated, the cell will have to generalize a lot of data values in order to determine its final color. If there is a large disparity in AOD throughout the region, this will not be visible if the cell color is based on the average or median reading level.
In this short vignette we will explore the effects of different raster grid sizes on collections of spatial points. A good place to find areas is Texas on May 9th, 2019, where the sea just off the east coast is dense with data points:
tb <- goesaodc_createTibble(nc) data_tx <- dplyr::filter(tb, lon > bbox_tx["x", "min"], lon < bbox_tx["x", "max"], lat > bbox_tx["y", "min"], lat < bbox_tx["y", "max"]) sp <- goesaodc_createSpatialPoints(nc, bbox = bbox_tx, dqfLevel = 2) par(mar=c(0, 0, 2, 0)) plot(tx, main = "Texas AOD") goesaodc_plotSpatialPoints(sp, var = "AOD", cex = 0.15, add = TRUE) plot(tx, add = TRUE)
Here are three rasters with different resolutions. Notice how the lower resolution grids suggest there is data for areas where no AOD measurements have been taken at all:
par(mfrow=c(1, 3), mar=c(2, 1, 2, 1)) rstr <- goesaodc_createRaster(nc, bbox = bbox_tx, res = 0.1, dqfLevel = 2) raster::plot(rstr$AOD, main = "AOD Res. 0.1", col = pal_aod(100), legend = FALSE) plot(tx, add = TRUE) rstr <- goesaodc_createRaster(nc, bbox = bbox_tx, res = 0.3, dqfLevel = 2) raster::plot(rstr$AOD, main = "AOD Res. 0.3", col = pal_aod(100), legend = FALSE) plot(tx, add = TRUE) rstr <- goesaodc_createRaster(nc, bbox = bbox_tx, res = 1.0, dqfLevel = 2) raster::plot(rstr$AOD, main = "AOD Res. 1.0", col = pal_aod(100), legend = FALSE) plot(tx, add = TRUE)
As the tiles increase in size they become more and more inaccurate in summarizing the true aerosol depth over the area they cover. Let's zoom to the dense cloud of points off the southeastern coast near Corpus Christi so we can better undsertand how data points influence individual cells. Each measurement location is marked with a black dot:
# Subset for the east coast of Texas where there is a lot of dense data bb_tx_dense <- c(-97, -96, 27, 28) tb_tx_dense <- dplyr::filter(tb, lon > bb_tx_dense[1], lon < bb_tx_dense[2], lat > bb_tx_dense[3], lat < bb_tx_dense[4]) rstr_mean <- goesaodc_createRaster(nc, bbox = bb_tx_dense, res = 0.1, fun = mean, dqfLevel = 2) # Plot the raster colored with the AOD palette and draw all sample point locations raster::plot(rstr_mean$AOD, main = "AOD", col = pal_aod(25), xlim = bb_tx_dense[1:2], ylim = bb_tx_dense[3:4]) plot(tx, add = TRUE) points(tb_tx_dense$lon, tb_tx_dense$lat, pch = 16, cex = 0.4)
Here we can see how the Advanced Baseline Imager (ABI) aboard the GOES-16 satellite takes its sample measurements in a non-orthogonal grid. Even in this uniform arrangement though there are still holes and gaps left in the readings. The cell colors in this plot are determined by the mean AOD value of all the points within the area they cover. In the dense spaces of the grid this may end up averaging around twenty points while the more more sparse regions to the upper-left may only use five or fewer values.
Drawing the sample points with the same AOD palette shows how the raster generalizes the individual readings within each cell. Let's compare how different our results are when using either the mean or median:
# Assign each reading its own color based on the AOD palette tb_tx_dense$col <- pal_aod(25)[as.numeric(cut(tb_tx_dense$AOD, breaks = 25))] # Generate the median raster rstr_median <- goesaodc_createRaster(nc, bbox = bb_tx_dense, res = 0.1, fun = median, dqfLevel = 2) par(mfrow = c(1, 2), mar = c(2, 4.5, 2, 4.5)) # Plot the mean raster raster::plot(rstr_mean$AOD, main = "Mean AOD", col = pal_aod(25), xlim = bb_tx_dense[1:2], ylim = bb_tx_dense[3:4]) plot(tx, add = TRUE) points(tb_tx_dense$lon, tb_tx_dense$lat, pch = 16, cex = 0.5, col = tb_tx_dense$col) # Plot the median raster raster::plot(rstr_median$AOD, main = "Median AOD", col = pal_aod(25), xlim = bb_tx_dense[1:2], ylim = bb_tx_dense[3:4]) plot(tx, add = TRUE) points(tb_tx_dense$lon, tb_tx_dense$lat, pch = 16, cex = 0.5, col = tb_tx_dense$col)
The two raster images turn out to be pretty similar. One could argue though that the median gives a better representation since the cells of that map blend better with the colors of all their points for some regions.
It's possible to use different summarizing statistics besides the mean and median to color cells. Rasters can be generated based on the minimum data value in a cell, the maximum value, the first or last, the total number of values, the standard deviation, and several other statistics. Let's try out the count and standard deviation options:
rstr_count <- goesaodc_createRaster(nc, bbox = bb_tx_dense, res = 0.1, fun = "count", dqfLevel = 2) rstr_sd <- goesaodc_createRaster(nc, bbox = bb_tx_dense, res = 0.1, fun = sd, dqfLevel = 2) par(mfrow = c(1, 2), mar = c(2, 4.5, 2, 4.5)) raster::plot(rstr_count$AOD, main = "Point Count", col = pal_count(25), xlim = bb_tx_dense[1:2], ylim = bb_tx_dense[3:4]) plot(tx, add = TRUE) points(tb_tx_dense$lon, tb_tx_dense$lat, pch = 16, cex = 0.3) raster::plot(rstr_sd$AOD, main = "AOD Standard Deviation", col = pal_sd(25), xlim = bb_tx_dense[1:2], ylim = bb_tx_dense[3:4]) plot(tx, add = TRUE, main = "AOD Standard Deviation") points(tb_tx_dense$lon, tb_tx_dense$lat, pch = 16, cex = 0.3)
Due to the irregularity of the data point layout, the count plot displays a sort of ripple effect. Meanwhile, the standard deviation plot is useful to determine areas where we may encounter the problem from earlier where a cell has to generalize values that differ greatly.
Now that we know what is going on in a small area, let's zoom back out to the state overview and see how various resolutions influence the apparent point counts and standard deviations of Texas:
par(mfrow=c(1, 3), mar=c(2, 1, 2, 1)) rstr <- goesaodc_createRaster(nc, bbox = bbox_tx, res = 0.1, fun = "count", dqfLevel = 2) raster::plot(rstr$AOD, main = "Point Count Res: 0.1", col = pal_count(100), legend = FALSE) plot(tx, add = TRUE) rstr <- goesaodc_createRaster(nc, bbox = bbox_tx, res = 0.3, fun = "count", dqfLevel = 2) raster::plot(rstr$AOD, main = "Point Count Res: 0.3", col = pal_count(100), legend = FALSE) plot(tx, add = TRUE) rstr <- goesaodc_createRaster(nc, bbox = bbox_tx, res = 1.0, fun = "count", dqfLevel = 2) raster::plot(rstr$AOD, main = "Point Count Res: 1.0", col = pal_count(100), legend = FALSE) plot(tx, add = TRUE)
par(mfrow=c(1, 3), mar=c(2, 1, 2, 1)) rstr <- goesaodc_createRaster(nc, bbox = bbox_tx, res = 0.1, fun = sd, dqfLevel = 2) raster::plot(rstr$AOD, main = "AOD SD Res: 0.1", col = pal_sd(100), legend = FALSE) plot(tx, add = TRUE) rstr <- goesaodc_createRaster(nc, bbox = bbox_tx, res = 0.3, fun = sd, dqfLevel = 2) raster::plot(rstr$AOD, main = "AOD SD Res: 0.3", col = pal_sd(100), legend = FALSE) plot(tx, add = TRUE) rstr <- goesaodc_createRaster(nc, bbox = bbox_tx, res = 1.0, fun = sd, dqfLevel = 2) raster::plot(rstr$AOD, main = "AOD SD Res: 1.0", col = pal_sd(100), legend = FALSE) plot(tx, add = TRUE)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.