Unsupervised Wildfire Detection Using GOES-16 Imagery

This notebook provides a Python pipeline of detecting wildfire from GOES-16 imagery using a robust principal component analysis (RPCA) based algorithm. The detail of this algorithm is described in Xu, Kong and Asgharzadeh (2021, link) which was presented in the 2021 IEEE International Symposium on Geoscience and Remote Sensing (IGARSS).

Table of Contents

Introduction¶

The Geostationary Operational Environmental Satellites (GOES)-16 [1], launched in 2017, is one of the newest weather satellites op-erated by NASA and the National Oceanic and AtmosphericAdministration (NOAA) to provide high frequency Earth imagery and atmospheric measurements. It scans the continen-tal U.S. (CONUS) every 5 minutes through 16 spectral bands covering visible and infrared wavelengths. The 2 km spatialresolution and high temporal resolution of GOES-16 imagery makes it suitable for detecting wildfires online.

Data¶

Real-time feed and historical archive of GOES-16 Imagery are publicly accessible through Amazon S3. The goes2go package provides function to download and visualize GOES imagery.

As an example, we can download the latest image and load it into an xarray.

The raw data are hyperspectral images across 16 bands. We can project them to a natural color RGB for visualization. We can also overlay the map with various features provided by cartopy. Here we specify coastlines and state boundaries.

For a full list of features, see http://www.naturalearthdata.com/downloads/

Lat/Lon Projection¶

The coordinates variables of GOES images are x an y which represent horizontal and vertical scan angles in radians.

Source: link

They can be used to reproject the images to Lat/Lon coordinates using the code below.

Now, Lat/Lon coordinates of (x,y) location can be obtained by

lat,lon = lats[x,y],lons[x,y]

We can also store the projected coordinates for future use.

gll = xr.Dataset(dict(latitude=(('y', 'x'), lats), 
                     longitude=(('y', 'x'), lons)))

gll.latitude.attrs = dict(standard_name='latitude', long_name='latitude', units='degrees_north')
gll.longitude.attrs = dict(standard_name='longitude', long_name='longitude', units='degrees_east')

gll.coords['latitude'] = gll.latitude
gll.coords['longitude'] = gll.longitude

gll = gll.where(~np.isinf(gll.latitude))
gll = gll.where(~np.isinf(gll.longitude))

ds.attrs['name'] = f'lat/lon grid for GOES-16 CONUS'

ds = ds.where(~np.isinf(ds.latitude))
ds = ds.where(~np.isinf(ds.longitude))

name = f'./latlon_GOES-16_CONUS.nc'
ds.to_netcdf(name)
ds

Kincade Fire¶

The Kincade fire was a wildfire that burned in Sonoma County, California in the United States. The fire started northeast of Geyserville in The Geysers on 9:24 p.m. on October 23, 2019 and subsequently burned 77,758 acres until the fire was fully contained on November 6, 2019. (source: wiki)

Fire Perimeter¶

Fire perimeter are GIS data that map the approximated location of fire. They serve as groundtruth for validating detection results. The National Interagency Fire Center provides fire perimeters for most wildfires that took place in 2019.

The variable sf is a list of shapeRecords, each shapeRecord stores the points (vertices) that form the fire perimeter.

We can convert the fire perimeter to a Polygon which is a spatial geometry. We can also visualize the Polygon on a map using the folium package.

We can check if a coordinate, represented by a Point object, falls within the Polygon.

Method¶

To detect Kincade Fire, we use GOES-16 ABI Radiance images from October 28.

Our method only requires 3.9 $\mathrm{\mu m}$ (IR shortwave window) and 12.3 $\mathrm{\mu m}$ ("Dirty" longwave window) band.

Each pixel in GOES-16 image represents a 2 km $\times$ 2 km region. For our analysis, we crop a 100 km $\times$ 100 km region in Northern California.

Radiance to Brightness Temperature Conversion¶

We convert radiance data to brightness temperature (BT) data using the function $$\text{BT}=\frac{\text{fk}_2}{\text{bc}_2\left[\log(\frac{\text{fk}_1}{L_{\lambda}}+1)-\text{bc}_1\right]}$$ where $L_{\lambda}$ is the radiance value and $\text{fk}_1,\text{fk}_2,\text{bc}_1,\text{bc}_2$ are wavelength dependent coefficients. Table of coefficients for different bands can be found in the GOES-R document.

Our detection statistic is the brightness temperature difference (BTD).

The BTD image can be seen as superposition of a spatially correlated background matrix and a sparse outlier matrix, where the outliers are potential fire pixels.

Sparse and Low-rank Matrix Decomposition¶

Robust Principal Component Analysis (RPCA) [2] decomposes an observation matrix into low-rank and sparse components through a two-objective optimization problem, $$\min_{\mathbf{L},\mathbf{S}}\ \lVert\mathbf{L}\rVert_*+\lambda\lVert\mathbf{S}\rVert_1\ s.t. \ \mathbf{Y}=\mathbf{L}+\mathbf{S}\label{eq2}$$ where $\lVert\mathbf{L}\rVert_*$ is the sum of eigenvalues $\sum_i\sigma_i(\mathbf{L})$, and $\lVert\mathbf{S}\rVert_1$ is the $\ell_1$-norm. The penalty coefficient $\lambda$ controls the sparsity of $\mathbf{S}$. It has been proved that under surprisingly broad conditions, the choice of penalty coefficient $\lambda=\frac{1}{\sqrt{\max(m,n)}}$ can exactly recover the low-rank and sparse components [2].

We use a python implementation of RPCA (source: https://github.com/dganguli/robust-pca).

Removing Cloud Pixels¶

The sparse matrix may contain cloud pixels which are false positives. Water vapor absorbs atmospheric energy at 12.3 $\mu$m, therefore clouds are significantly cooler than the background in the 12.3 $\mu$m band. We propose to identify a pixel in $\mathbf{S}$ as a cloud pixel if its radiance at the 12.3 $\mu$m band, $R_{12.3}$, is significantly low.

The T-point algorithm [3] is an automatic image thresholding method based on the pixel density function. It assumes that the pixel density can be approximated by a unimodal histogram with a heavy upper tail. The histogram can be decomposed into three parts: a steep rising slope, a steep descending slope, and a slow descending slope. It then fits a piece-wise linear regression to the steep descending and the low descending slope. The optimal cutoff value is set equivalent to the knot that minimizes the fitting error, such as the $\ell_d$-norm, $d\in\{1,2\}$, of the regression. Let $\mathbf{r}$ denote a vector that stores the sorted negative pixel $R_{12.3}$ of an $m$ by $n$ image in ascending order. If cloud pixels contribute to the lower tail of the $R_{12.3}$ density, then they will contribute to the higher tail of the $-R_{12.3}$ density. Let $\mathbf{r}_{(i)}, \ i\in[1,mn]$ be the $-R_{12.3}$ of the $i$-th entry. If we assume the peak of $-R_{12.3}$ histogram has been identified at index $p$, then the optimal cutoff index $t\in[p+1,mn-1]$ is found by solving the optimization problem $$\begin{equation} t = \text{arg min}\left\{ \sum_{i=p}^{t}\left\lvert \mathbf{r}_{(i)}-\hat{\mathbf{r}}_{(i)}\right\rvert^d+\sum_{i=t}^{mn}\left\lvert \mathbf{r}_{(i)}-\hat{\mathbf{r}}_{(i)}\right\rvert^d\right\} \label{eq3}\tag{1} \end{equation}$$ where $\hat{\mathbf{r}}$ are the estimated $-R_{12.3}$ by the fitted regression lines. Following Equation $\eqref{eq3}$, we will identify a pixel as wildfire only if its $R_{12.3}$ is above $-\mathbf{r}_{(t)}$.

Removing Noise Pixels¶

The sparse matrix can also contain noise due to sun glints and RPCA algorithm itself. We apply a low-pass filter and rule-based threshold to attenuate high-frequency noise.

A threshold value between 5 and 8 generally performs well.

Result¶

We can overlay our detection on the map we created before.

Our method detected two wildfires, the left one is the Kincade Fire as confirmed by the fire perimeter, the right one is actually the Grizzly Island Fire which started on Oct 28th.

Comparison with GOES-16 Active Fire Product¶

GOES-16 has a built-in fire detection algorithm [4]. Its results can be extracted similar like other products.

The FDC algorithm assigns each pixel a mask code. Among them, 10-15 and 30-35 correspond to active fire (source: link).

We see that although FDC successfully identified the two wildfires, it is over-optimistic about the severity of Kincade Fire and yields significantly smaller coverage. In contrast, our method yields a more accurate characterization of the current status of the fire and can help stakeholders better estimate the risk and resources needed to suppress the fire.

Reference ¶

[1] GOES-R Calibration Working Group and GOES-R Series Program, (2017): NOAA GOES-R Series Advanced Baseline Imager (ABI) Level 1b Radiances. NOAA National Centers for Environmental Information. doi:10.7289/V5BV7DSR.

[2] Emmanuel J Candès, Xiaodong Li, Yi Ma, and John Wright, "Robust principal component analysis?," Journal of the ACM (JACM), vol. 58, no. 3, pp. 1–37, 2011.

[3] Nicolas Coudray, Jean-Luc Buessler, and Jean-Philippe Urban, "Robust threshold estimation for images with uni-modal histograms," Pattern Recognition Letters, vol. 31,no. 9, pp. 1010–1019, 2010.

[4] Christopher C. Schmidt, Jay Hoffman, Elaine Prins, and Scott Lindstrom (2013): GOES-R Advanced Baseline Imager (ABI) Algorithm Theoretical Basis Document For Fire / Hot Spot Characterization. NOAA NESDIS Center for Satellite Applications and Research.