Estimating Atmospheric Motion Winds from Satellite Image Data using Space-time Drift Models

by   Indranil Sahoo, et al.

Geostationary satellites collect high-resolution weather data comprising a series of images which can be used to estimate wind speed and direction at different altitudes. The Derived Motion Winds (DMW) Algorithm is commonly used to process these data and estimate atmospheric winds by tracking features in images taken by the GOES-R series of the NOAA geostationary meteorological satellites. However, the wind estimates from the DMW Algorithm are sparse and do not come with uncertainty measures. This motivates us to statistically model wind motions as a spatial process drifting in time. We propose a covariance function that depends on spatial and temporal lags and a drift parameter to capture the wind speed and wind direction. We estimate the parameters by local maximum likelihood. Our method allows us to compute standard errors of the estimates, enabling spatial smoothing of the estimates using a Gaussian kernel weighted by the inverses of the estimated variances. We conduct extensive simulation studies to determine the situations where our method should perform well. The proposed method is applied to the GOES-15 brightness temperature data over Colorado and reduces prediction error of brightness temperature compared to the DMW Algorithm.



There are no comments yet.


page 16

page 18

page 24


Sensitivity Analysis of Wind Energy Resources with Bayesian non-Gaussian and nonstationary Functional ANOVA

The transition from non-renewable to renewable energies represents a glo...

Scene Learning: Deep Convolutional Networks For Wind Power Prediction by Embedding Turbines into Grid Space

Wind power prediction is of vital importance in wind power utilization. ...

Gaussian process aided function comparison using noisy scattered data

This work proposes a new nonparametric method to compare the underlying ...

Measuring Wind Turbine Health Using Drifting Concepts

Time series processing is an essential aspect of wind turbine health mon...

Multi Scale Graph Wavenet for Wind Speed Forecasting

Geometric deep learning has gained tremendous attention in both academia...

Using neural networks to estimate parameters in spatial point process models

In this paper, I show how neural networks can be used to simultaneously ...

Log-Gaussian Cox Process Modeling of Large Spatial Lightning Data using Spectral and Laplace Approximations

Lightning is a destructive and highly visible product of severe storms, ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Winds are an important component of atmospheric circulation, and estimating local or regional winds is particularly important for making weather forecasts. They also influence vegetation as they affect factors that influence plant growth, such as seed dispersal rates, air transportation of pollen, and metabolism rates in plants. The study of local winds also permits evaluation of power produced by wind turbines (Brown et al., 1984; Castino et al., 1998), prediction of propagation of oil-spills (Kim et al., 2014) and the study of coastal erosion (Ahmad et al., 2015). Local winds are also capable of moving pollutants into an area. For example, Calima, which blows dust into the Canary Islands (WeatherOnline, 2018). Winds may also impact large-scale devastation such as forest fires. For example, the Santa Ana winds which blow into California after scorching summers (Berkowitz & Steckelberg, 2017). Thus it is important to map the strength and direction of local winds to help us prepare for natural calamities and facilitate preventive measures.

Observations of wind speed and direction at the ground level are collected at weather stations on land and by buoys or ships over oceans. Low orbit satellites such as Jason 3 and Sentinel 1 infer surface winds using geophysical inversion algorithms, based on peak backscattered power and the shape of radio signal waveforms (ESA, 2016). However, ground monitors are sparse in space, and the satellites cannot monitor winds continuously in space and time as they are in low earth orbit. Winds in the upper level of the atmosphere can be observed using weather balloons or aircraft measurements, but these observations are also very sparse in space and time.

There is a significant statistical literature on modeling winds from ground monitors (Priestley, 1981; Haslett & Raftery, 1989; Brillinger, 2001; Stein, 2005). For some applications, it is sufficient to look at the evolution of wind at a fixed location, and several methods have been proposed for this scenario (Brown et al., 1984; Tol, 1997; Ailliot, 2004; Monbet et al., 2007). The wind at different locations can be utilized to model spatiotemporal dependencies (Bennett, 1979; Bras & Rodriguez-Iturbe, 1985; Kyriakidis & Journel, 1999; De Luna & Genton, 2005). Also, Boukhanovsky et al. (2003), Malmberg et al. (2005) and Ailliot et al. (2006) proposed autoregressive space-time models to describe the evolution of winds. Stein (2005) proposed a spectral-in-time modeling approach to describe the space-time dependencies of the data. Fuentes et al. (2008) modeled a drift process using Bayesian analysis, where the drift parameter is modeled using splines. Modlin et al. (2012)

used circular conditional autoregressive models for wind direction and speed.

On the other hand, geostationary weather satellites provide data from the surface and the atmosphere with a very high temporal resolution. The resulting data comprise a series of images which essentially make them a ‘movie’. While the satellites do not directly measure wind, the image sequences are used to infer wind estimates by tracking movements of atmospheric tracers such as clouds or moisture features over time. Wind data obtained from satellite images play a major role in data assimilation. Numerical climate models perform better with accurate wind data, especially over the oceans, resulting in improved weather forecasts and warnings (Tomassini et al., 1999). For example, the European Centre for Medium-Range Weather Forecasts (ECMWF) has been incorporating atmospheric motion winds into their forecast models operationally since the 1980s. This has dramatically improved the model’s ability to forecast the track of tropical cyclones and has also increased the model’s ability to predict wave heights and storm surges (Tomassini et al., 1999).

The Derived Motion Winds (DMW) Algorithm (Daniels et al., 2010) is a standard algorithm for estimating motion winds from satellite images. The DMW Algorithm takes as input brightness temperature images from the National Oceanic and Atmospheric Administration (NOAA) geostationary meteorological satellites and gives estimated wind fields as outputs. The algorithm tracks a suitable target across the input images and assigns a motion wind to the middle time point (see Section 2 for details). The image at the middle time point must satisfy a set of criteria to qualify as a suitable target scene. As a result, the DMW estimates are often missing. The DMW algorithm also does not produce a measure of uncertainty.

Geostatistical methods are capable of overcoming these shortcomings. This motivates us to model satellite image data using a spatial process drifting in time. At the heart of this statistical model lies the idea of incorporating the motion vector parameters in the process covariance.

Stein et al. (2013) uses this idea to fit a space-time model. We borrow the idea of Nested Tracking from Daniels et al. (2010) by considering data buffers and estimate local wind vectors using maximum likelihood estimates over sliding windows over space. Local estimation of covariance parameters using moving windows has been studied in Haas (1990, 1995). One major advantage of our approach over the DMW algorithm is that it allows us to quantify uncertainties associated with the estimates. The estimated wind fields are smoothed using weighted Gaussian kernels, the kernel being scaled by these estimated inverse variances. This not only enables us to make smooth maps of wind over specific regions but also improves forecast of brightness temperature fields, which is evidence that the wind fields are estimated more accurately.

2 Derived Motion Winds Algorithm

The Derived Motion Winds (DMW) Algorithm estimates atmospheric motion winds from images taken by geostationary satellites. For a cloudy region, the imager records brightness temperature (see Figure 1), which measures the radiance (in Kelvin) of microwave radiation traveling upward from the top of the atmosphere to the satellite. For clear sky portions, the satellite records images of suitable indicators of atmospheric moisture content, such as specific humidity. Daniels et al. (2010) provides a description of and the physical basis for the estimation of atmospheric winds from the images taken by geostationary satellites.

The DMW algorithm involves creating a data buffer

, which is a data structure holding 2-dimensional arrays of the response variable for 3 consecutive image times. The middle portion of the buffer is divided into smaller ‘target scenes’, and each scene is analyzed to locate and select a set of suitable targets in the middle image.

Daniels et al. (2010) gives a description of the Nested Tracking Algorithm which involves nesting smaller target scenes (usually of size ) within a large target scene of size pixels and getting every possible local motion vectors derived from each possible smaller box within a large target scene. The displacement vector between time points and is computed by minimizing the Sum of Squared Differences (SSD) criterion as

where denotes the brightness temperature within the smaller box at pixel location and time point , is the indices of pixels in the scene centered at , and is a two-dimensional vector denoting the displacement. The sum is considered over two dimensions and the optimization over is done only over integers so that corresponds to an observed pixel. In practice, the search region is substantially larger than the size of the smaller target scene, so the above summation is carried out for all target box positions within the search region. The mean displacement vector is computed as

and is assigned as the DMW estimate at location time point

in the buffer. Once every possible local motion vectors within the buffer are calculated, a density-based cluster analysis algorithm, DBSCAN

(Ester et al., 1996) is used to identify the largest cluster representing the dominant motion. The final DMW estimate for the buffer is the average of the vectors belonging to the largest cluster.

The size of the target scene depends on the spatial and temporal resolution of the imagery and the scale of the intended feature to be tracked. Daniels et al. (2010) suggests that the temporal resolution of the images should at most be 15 minutes in order to account for the short lifespan and rapid disintegration of clouds over land. The DMWA does not offer wind estimates at every as the data in the middle image has to satisfy a set of criteria to qualify as a suitable target scene. Also, quantifying uncertainties using the SSD criterion is hard because each estimated vector uses a different subset of the data, so likelihood ratio tests are not applicable. Finally, the vector estimates generated for each target scene can at most be half-integers.

3 Model-based wind estimation

3.1 Space-time drift models

The proposed approach uses spatiotemporal covariance functions to track the wind. This requires us to consider asymmetric spatiotemporal covariance functions. The space-time process has asymmetric covariance if


In most regions, winds flow in a consistent direction, and so changes in temperature or precipitation rate at one location tend to precede similar changes down wind. For instance, if and winds flow consistently from to , then we expect

We incorporate space-time asymmetries via a drift parameter. Suppose that is a stationary, space-time symmetric process with covariance function

and let

Then the covariance function of is


which is space-time asymmetric and stationary. The parameter can be interpreted as the drift of the process over time, which we use to estimate winds.

3.2 Local estimation of the drift parameter

Assume the brightness temperature data have been standardized to have mean zero and variance one at each location (as described in the Appendix) and denote the standardized data as . The mean and variance carry no information about the drift, and this step simplifies estimation of the covariance parameters. We do not specify the global covariance function for . Instead, we specify its local covariance with drift parameter and estimate the motion winds locally.

We define a target scene as a square array of pixels

To estimate , we assume that the process is locally stationary in and that winds are smooth enough to be assumed constant in the scene, that is, . In other words, we approximate the local covariance function as

where denotes the spatial lag, denotes the temporal lag and . In particular, we assume


Here, and denote respectively the spatial and temporal range parameters. Let be the four correlation parameters to be estimated.

We use maximum likelihood estimation within to estimate . That is, if denotes the standardized data vector in the scene and denote the corresponding space-time covariance matrix with elements defined by (3), then the log-likelihood for given is


The estimates obtained are associated with the location and time point at which the target scene was centered, denoted . We also estimate the variances associated with the estimated wind vectors by computing the observed inverse Hessian matrix at the maximum likelihood estimate. We imitate the Nested Tracking approach and slide the scene window across space and time, estimating wind vectors locally in space and time using the same optimization routine.

3.3 Smoothing the local estimates

After obtaining the local estimates of for all , we smooth these initial estimates to stabilize them by borrowing strength across space. The two components of the wind vectors are smoothed separately. The kernel smoothing weights are taken to be proportional to the ratio of a spatial Gaussian kernel and the variance of the initial estimate. Full details are given in the Appendix. The bandwidth is chosen based on cross validation.

The size of the target scene is an important tuning parameter in the study. While implementing the method on real data sets, the window size should be chosen such that the wind motion is roughly constant in the scene, and the feature being tracked in time is prominent and does not move out of frame. In Section 4, we perform a simulation study analyzing the effect of window size on the performance of our method. In the real-data analysis in Section 5, the optimal window size is chosen using cross validation.

4 Simulation results

In this section, we conduct a simulation study to determine the conditions under which the space-time drift model (STDM) performs well for estimating wind motion vectors. For this purpose, we repeatedly simulate datasets within one particular target scene (as opposed to scanning across a spatial domain) and we also implement a version of the DMW algorithm and compare its performance with the STDM. To compare the two methods, accuracy for simulated dataset is measured by Vector Difference (Daniels et al., 2010) between the true () and estimated () wind vectors

and we report the Mean Vector Difference over datasets

and the standard deviation in Table


First we consider target scenes of size generated independently from the STDM under various parameter settings. The true spatial range parameter is chosen to be either 1, 2, 4 or 8 and the true temporal range parameter is chosen to be either 1, 2, 3 or 4. We also take two different values of the reference wind vector, namely and , which signify respectively slow and fast wind vectors. All four parameters are updated simultaneously during optimization. Table 1 shows the performance of the two methods for the two wind vectors based on simulations.

MVD of STDM for
1 2 3 4
1 0.415 (0.37) 0.172 (0.09) 0.136 (0.07) 0.118 (0.05)
2 1.225 (1.21) 0.304 (0.19) 0.162 (0.09) 0.149 (0.08)
4 3.010 (2.56) 1.008 (0.72) 0.398 (0.25) 0.268 (0.14)
8 3.441 (3.48) 2.830 (2.06) 1.346 (1.00) 0.831 (0.69)
MVD of DMWA for
1 2 3 4
1 1.965 (1.17) 0.771 (0.91) 0.162 (0.50) 0.084 (0.42)
2 2.230 (1.30) 1.343 (1.08) 0.727 (0.82) 0.409 (0.69)
4 2.532 (1.32) 1.888 (1.18) 1.856 (1.17) 1.101 (0.86)
8 2.891 (1.54) 2.538 (1.33) 2.064 (1.09) 1.961 (1.19)
Coverage of STDM for
1 2 3 4
1 84 91 84 85
2 79 93 90 88
4 59 81 91 93
8 65 64 79 95
MVD of STDM for
1 2 3 4
1 1.124 (1.09) 0.868 (1.50) 0.743 (1.57) 0.530 (1.06)
2 1.916 (1.69) 0.777 (1.24) 0.299 (0.42) 0.230 (0.41)
4 2.820 (1.98) 1.489 (1.36) 0.676 (0.70) 0.392 (0.37)
8 3.401 (2.84) 3.125 (1.96) 2.071 (1.80) 1.129 (1.03)
MVD of DMWA for
1 2 3 4
1 5.823 (1.62) 5.836 (1.43) 5.866 (1.59) 6.024 (1.61)
2 5.638 (1.74) 5.294 (1.75) 5.408 (1.80) 5.362 (1.64)
4 5.548 (1.42) 4.995 (1.65) 4.934 (1.85) 4.750 (1.70)
8 5.864 (1.62) 5.320 (1.72) 5.233 (1.75) 4.603 (1.70)
Coverage of STDM for
1 2 3 4
1 73 77 78 78
2 66 89 86 88
4 63 72 86 87
8 66 64 76 84
Table 1: Comparing Mean Vector Difference (SD) for Space-time Drift Model (STDM, first row) and Derived Motion Winds algorithm (DMWA, second row) for true wind vectors (left panel) and (right panel) based on data windows of size ; and denote the true spatial and temporal range respectively. The third row gives the 95% coverage of STDM for the two reference vectors, averaged over the coverage for the u- and v- components.

STDM does a better job in estimating moderately small and large wind vectors for the window size compared to DMWA. Mean vector distance is the smallest when the true spatial range is small and the true temporal range is large. This is intuitive because a small spatial range makes it easier to identify a feature in the target scene, and a large temporal range means that the features dissipate slowly over time. STDM also performs better for the smaller wind vector because when the wind vector is large compared to the window size, the feature tracked in time could potentially move out of the frame, resulting in incorrect wind estimates.

The performance of the DMW algorithm also improves as the temporal range increases. However, these simulation results bring forth a major flaw in the DMW algorithm. The estimated motion winds from the DMW algorithm are at most half integers and they are limited to the size of the larger search window. That is, while estimating the motion vectors, the smaller central target scene () can only move up to 4 pixels in all directions while it is being tracked back and forward in time. As a result, it performs poorly for the larger wind motion vector, which had a v-component of 5 pixel units. This can once again be attributed to the window size relative to the magnitude of the wind vector.

To examine the effect of window size, we perform the simulations again under the same covariance parameter settings as described earlier but with window sizes of and respectively. Tables 1, 2 and 3 give us a clear picture of the effect of window size on the estimation. Target scenes of size are not adequate to contain the features being tracked in frame for a wind vector of which is reflected in the high MVD and SD values for both STDM and DMWA (see Table 2). The estimation improves as we increase window size to (see Table 1) and then to (see Table 3). Tables 1, 2 and 3 also show that the estimation of the wind vector improves with the increase in grid size.

Another way to assess the effect of window size on STDM is to look at coverage probabilities for the u- and v- components of the wind vector. Since our model allows uncertainty quantification through estimated variances of the estimates, we can form confidence intervals for the parameters of interest. This exploits the property that the maximum likelihood estimates are asymptotically normal under increasing domain. Once again we consider

to be either or and look at 95% coverage probabilities of the two components for window sizes , and respectively. The other covariance parameters are the same as before. For each scenario, the coverage probability is computed based on 100 replications. The coverage for the three different scenarios, averaged over the u- and v-components have been shown in Tables 1 - 3. Detailed version of these results are shown in the Supplementary Materials.

MVD of STDM for
1 2 3 4
1 1.196 (0.35) 0.311 (0.22) 0.243 (0.24) 0.201 (0.20)
2 1.658 (1.13) 0.600 (0.52) 0.353 (0.33) 0.245 (0.12)
4 2.841 (1.75) 1.599 (1.24) 0.861 (0.61) 0.456 (0.31)
8 3.253 (2.25) 3.281 (1.87) 2.212 (1.71) 1.432 (1.02)
MVD of DMWA for
1 2 3 4
1 1.983 (1.05) 1.209 (0.85) 1.039 (1.12) 0.854 (1.09)
2 2.045 (1.02) 1.527 (0.94) 1.283 (1.01) 1.072 (0.98)
4 2.342 (0.94) 2.158 (0.94) 1.620 (0.92) 1.511 (0.91)
8 2.452 (1.07) 2.175 (0.95) 2.003 (1.10) 1.776 (0.89)
Coverage of STDM for
1 2 3 4
1 66 90 93 96
2 72 88 91 95
4 74 76 91 97
8 61 85 90 87
MVD of STDM for
1 2 3 4
1 2.310 (1.60) 2.034 (1.66) 1.734 (1.68) 1.787 (1.95)
2 2.335 (1.56) 2.032 (1.76) 1.392 (1.47) 1.030 (1.21)
4 3.440 (2.34) 2.393 (1.90) 1.996 (1.70) 1.638 (1.99)
8 3.611 (2.66) 3.135 (1.92) 2.913 (1.88) 2.483 (2.39)
MVD of DMWA for
1 2 3 4
1 5.953 (0.91) 5.924 (0.96) 5.953 (1.10) 6.012 (0.98)
2 5.710 (1.09) 5.781 (1.18) 5.899 (1.09) 5.739 (1.08)
4 5.665 (1.09) 5.577 (1.16) 5.540 (1.21) 5.355 (1.13)
8 5.952 (1.13) 5.678 (1.23) 5.413 (1.12) 5.461 (1.18)
Coverage of STDM for
1 2 3 4
1 51 59 55 62
2 43 49 64 68
4 55 59 54 59
8 73 61 57 60
Table 2: Comparing Mean Vector Difference (SD) for Space-time Drift Model (STDM, first row) and Derived Motion Winds algorithm (DMWA, second row) for true wind vectors (left panel) and (right panel) based on data windows of size ; and denote the spatial and temporal range respectively. The third row gives the 95% coverage of STDM for the two reference vectors, averaged over the coverage for the u- and v- components.

The coverage probabilities reiterate the conditions that we had arrived at from the MVD values in Tables 1 - 3. In particular, our model has better coverage when the spatial range is small and the temporal range is large. Once again, we see the importance of the window size. For instance, the coverage probability of our model is very small (around 60%) when we use a window size of to estimate a large wind of , showing that the window is not adequate to track the feature across time. Coverage increases as the grid size increases and achieves the nominal level for sufficiently large window size (Table 3).

In this simulation, a larger target scene always produces the most accurate estimates of the wind components. However, this might not be the case in practice. In the simulation study, the reference winds have been chosen to be uniform across space and time. In practice, this assumption might not be valid and choosing an arbitrarily large target scene might lead to inaccurate estimates of the wind vectors. Therefore, selecting the window size requires balancing a tradeoff between a window size that is large enough to capture targets moving through the scene yet small enough to satisfy the assumption that within the window the process is stationary with a constant drift. While implementing the methods on real data sets, we choose the optimal window size using cross validation. This should work well as wind fields are mostly smooth over a small region and time frame.

MVD of STDM for
1 2 3 4
1 0.274 (0.32) 0.125 (0.06) 0.097 (0.05) 0.073 (0.04)
2 0.803 (0.71) 0.197 (0.11) 0.122 (0.07) 0.103 (0.05)
4 2.443 (1.64) 0.536 (0.40) 0.251 (0.14) 0.183 (0.08)
8 3.268 (3.60) 2.030 (2.05) 0.875 (0.59) 0.496 (0.35)
MVD of DMWA for
1 2 3 4
1 2.998 (1.55) 1.501 (1.56) 0.699 (1.33) 0.318 (0.84)
2 2.776 (1.42) 2.028 (1.40) 1.365 (1.57) 0.815 (1.18)
4 3.213 (1.69) 2.661 (1.67) 2.236 (1.67) 1.644 (1.47)
8 3.579 (1.70) 2.921 (1.39) 2.700 (1.53) 2.518 (1.34)
Coverage of STDM for
1 2 3 4
1 91 95 93 93
2 86 90 92 93
4 71 91 93 93
8 67 78 87 88
MVD of STDM for
1 2 3 4
1 0.581 (1.16) 0.403 (1.21) 0.355 (0.88) 0.291 (0.99)
2 1.073 (1.14) 0.244 (0.13) 0.142 (0.07) 0.107 (0.06)
4 2.781 (2.21) 0.690 (0.65) 0.307 (0.18) 0.211 (0.11)
8 3.175 (3.46) 2.721 (2.55) 1.123 (0.97) 0.665 (0.53)
MVD of DMWA for
1 2 3 4
1 5.100 (2.54) 2.830 (2.46) 0.688 (1.54) 0.098 (0.62)
2 5.306 (2.43) 3.417 (2.51) 1.819 (2.16) 1.028 (1.77)
4 5.724 (2.47) 3.951 (2.42) 3.378 (2.37) 2.670 (2.32)
8 6.109 (2.26) 5.15 (2.44) 4.393 (2.33) 4.061 (2.56)
Coverage of STDM for
1 2 3 4
1 85 90 84 85
2 84 92 93 93
4 71 84 92 91
8 62 67 78 88
Table 3: Comparing Mean Vector Difference (SD) for Space-time Drift Model (STDM, first row) and Derived Motion Winds algorithm (DMWA, second row) for true wind vectors (left panel) and (right panel) based on data windows of size ; and denote the spatial and temporal range respectively. The third row gives the 95% coverage of STDM for the two reference vectors, averaged over the coverage for the u- and v- components.

5 Application to GOES-15 data

5.1 GOES-15 data description

The Geostationary Operational Environmental Satellite (GOES), operated by NOAA provides continual measurements of the atmosphere and surface variables, which help facilitate meteorological research including weather forecasting and severe storm tracking. Since the launch of GOES-8 in 1994, the GOES instruments have monitored atmospheric phenomena and provided a continuous stream of environmental data. The dataset used in this project is from the GOES-15 satellite. Launched in March 2010, GOES-15 is positioned at the GOES-West location of 135°W longitudes over the Pacific Ocean. The dataset, as described in Knapp & Wilkins (2018), is a gridded satellite Contiguous US domain data which are geostationary data remapped to equal angle projection with an 0.04°(4 km) latitudinal resolution and 15 minutes temporal resolution. The dataset includes infrared channel data in terms of pixel-wise brightness temperature for the reflective bands (channels 1 - 6 with approximate central wavelengths 0.47, 0.64, 0.865, 1.378, 1.61, 2.25 microns respectively). The reflective bands support among other ground and atmospheric indicators, the characterization of clouds. Gridded GOES-15 data can be obtained from NOAA One-Stop at

We analyze Channel 4 brightness temperature data (recorded in Kelvin scale) for 10 consecutive days starting January 1, 2015 at a temporal resolution of 15 minutes (960 total images), covering a region in Colorado ( N to N latitudes and W to W longitudes). For our analysis, we have focused on the region of Northeast Colorado, comprising 6,160 pixels, where we can see some clear cloud movements from south-west to north-east. Figure 1 shows the data at three consecutive time points on the third day. Lower values of brightness temperature indicates cloud cover, whereas high brightness temperature values suggest clear skies over the region.

Figure 1: Brightness temperature (Kelvin) maps over Colorado on January 3, 2015 at 00:52 am, 01:07 am and 01:22 am respectively.

5.2 Estimation using GOES-15 data

The first step involves standardizing the brightness temperature data using the pixel-wise sample mean and standard deviation over time. This factors out the effect of low cloud cover over the region which can affect the local estimation of wind vectors. To smooth out the standard deviation map, we use a Gaussian kernel with smoothing parameter pixels (i.e., 8 km).

We use cross validation to determine the appropriate size of target scenes. With no direct measurement of the wind, we compare window sizes indirectly based on their predictions of brightness temperature at the fourth time point using the wind estimate based on the first three time points and use this optimal window size to estimate winds for the subsequent time points. Table 4 gives us the Mean Squared Prediction Error (MSPE) (discussed in Section 5.3) for a few window sizes and the corresponding computation time to estimate the wind vectors at all 6,160 pixels using 4 consecutive time steps. Based on these results, we use window size of pixels (i.e., 100 km 100 km regions) and estimate the wind vectors locally at each spatial location and at each time using maximum likelihood estimates using the data in the window. We also estimate the variances associated with the estimates by computing the inverse Hessian matrix at the MLE. Figure 2 shows the uncertainty associated with the estimated wind components as given by the estimated standard deviations in the log scale, for three consecutive time points.

window size MSPE Computation time (in hours)
11 0.409 1.083
15 0.276 3.371
21 0.225 12.76
25 0.208 24.85
35 0.214 81.00
Table 4: Mean Squared Prediction Error for cross validation based on the first three time points and the corresponding computation time (in hours)
Figure 2: Estimated standard deviations (in log scale) obtained using the Space-Time Drift Model (STDM), corresponding to estimated u- (top row) and v- (bottom row) components at three consecutive time points. The columns (from left to right) represent 3 consecutive time points, t = 2, 3 and 4 respectively.

The estimated standard deviations in Figure 2 are high for some locations, perhaps due to convergence problems, which makes the wind estimates rough. To account for this, we smooth each component of the estimated wind field using weighted Gaussian kernels, the weights being scaled to the inverse variances of the estimates (details are given in the Appendix). The smoothing parameter has been chosen based on cross validation. Figure 3 shows the raw and smoothed wind estimates obtained from the proposed STDM. From Figures 2 and 3, it can be seen that the rough wind estimates correspond to regions where the estimated standard deviation is high. Figure 3 also shows the smoothed estimates. For fair comparison, the DMW estimates are also smoothed using a simple Gaussian kernel smoother with its optimal smoothing parameter chosen using cross validation. The raw and smoothed estimates of the wind fields obtained from DMWA are shown in Figure 4.

Figure 3: Raw (top row) and smoothed (bottom row) wind field estimates at three consecutive time points obtained using the Space-Time Drift Model (STDM). The columns (from left to right) represent 3 consecutive time points, t = 2, 3 and 4 respectively. The background is the county boundaries in northeast Colorado.
Figure 4: Raw (top row) and smoothed (bottom row) wind field estimates at three consecutive time points, obtained from the Derived Motion Winds Algorithm (DMWA). The columns (from left to right) represent 3 consecutive time points, t = 2, 3 and 4 respectively. The background is the county boundaries in northeast Colorado.

5.3 Comparison based on Mean Squared Prediction Error

For this dataset, reference wind fields are not available. Hence, we compare the two methods based on Mean Squared Prediction Error (MSPE) while predicting standardized brightness temperature fields. To predict at a spatial location at time , we consider , the standardized data in and estimated winds at time . This is because the estimated winds at time uses data from time . The predicted standardized temperature is calculated as the mean of the Gaussian conditional distribution of given under the stationary space-time drift model with estimated drift parameter . We also consider a naive approach of predicting the brightness temperature fields, where the the data at time is considered to be the predicted standardized brightness temperature fields at time . We call it the baseline prediction and the performances of the two methods are assessed relative to the baseline MSPE. Table 5 compares the raw and smoothed estimates of the wind components in terms of MSPE.

STDM 0.278 0.195 0.205 0.224
DMWA 0.400 0.293 0.298 0.265
Smoothed STDM 0.278 0.195 0.180 0.194
Smoothed DMWA 0.364 0.251 0.258 0.214
Baseline 1.275 1.180 0.928 0.812
Table 5: Comparing Mean Squared Prediction Error based on prediction using the raw and smoothed wind estimates from the Space-Time Drift Model (STDM) and the Derived Motion Winds Algorithm (DMWA) at different time points. in each case denotes the optimal smoothing parameter chosen using cross validation. All the methods have been compared against the baseline prediction.

From Table 5, we can conclude that the proposed model outperforms the DMW algorithm. We also see that smoothing the estimates does provide a better picture of the wind fields over Northeast Colorado. Figure 5 provides maps of predicted standardized brightness temperature using the estimated winds from STDM along with the lower and upper prediction maps for the brightness temperature fields, calculated using prediction variances. Figure 5 also shows the observed temperature fields for the three consecutive time points. Figure 5 shows that our model captures the main features in the temperature fields over Northeast Colorado. Our model also captures the wind movement over the region since we can see similar movements of features across the region as compared to the ones tracked along in the original images.

Figure 5: The top row shows predicted (standardized) brightness temperature fields for three consecutive time points, obtained using smoothed wind estimates from the Space-Time Drift Model (STDM). The second and third rows show respectively, the corresponding lower and upper prediction regions. The bottom row shows the observed standardized temperature fields at those time points.

6 Discussions and Conclusions

Wind is one of the most important atmospheric variables; it has large impact on local weather and hence, studying winds is essential. Wind data can be derived from high-resolution spatiotemporal data collected by geostationary satellites. These data are sequence of images over time and are used to derive atmospheric wind speed and direction. One algorithm that provides wind estimates is known as the Derived Motion Winds Algorithm. It takes as its input a sequence of images taken by the GOES-R series of the NOAA meteorological satellites and produces wind speed and direction. However, this algorithm does not quantify uncertainties. In this paper, we propose a spatiotemporal model to analyze satellite image data with the primary objective of estimating atmospheric wind speed and direction.

We have proposed local estimation of drift parameters. Developing a globally valid nonstationary space time drift model is an interesting problem that we have not pursued here. Following the basic idea of Nested Tracking, we propose a method to estimate the covariance parameters using maximum likelihood estimation. We smooth the raw estimated wind fields using a weighted Gaussian kernel, the weights being scaled by the inverse of the estimated variances of the estimates. Section 4 details an extensive simulation study that outlines conditions under which our model performs well. Based on our simulation study, we conclude that we have accurate wind estimates when the true spatial correlation range is small and the true temporal correlation range is high. The simulations also bring forth a major drawback of the DMWA. Due to the design of the DMW algorithm, the local DMWA estimates can only be half-integers and can only take values equal to the number of pixels the smaller target scene can move around in the larger search window. This brings us to perhaps the most important tuning parameter in the analysis, the target window size. We show that the window size is very crucial for both methods with large bias resulting from a large window and variance resulting form a small window.

We apply our method on brightness temperature data over Northeast Colorado, obtained from the GOES-15 satellite. While estimating winds, the window size has been chosen using cross validation. We provide estimated standard deviation maps, driving home the point that our method is capable of quantifying uncertainties associated with the estimation. We also provide smoothed maps of estimated wind fields over Northeast Colorado. We predict brightness temperature fields using our model and conclude, based on Mean Squared Prediction Error (MSPE) that the smoothed version of the wind estimates better represent the true wind condition. We also compare the performance of our proposed method and the DMWA. We have shown that our method outperforms the DMWA with respect to MSPE. We argue that smoothing the estimated wind fields give more reliable wind estimates. We also see that we capture the main features in the brightness temperature maps through our prediction, including the drift across the region.

One of the major challenges of our method is to apply it to data in real time. The main computational bottleneck is the time required for optimization over the covariance parameters. Because we use a local likelihood approach, the method is embarrassingly parallelizable across pixels and time points, and thus should scale well when used in production. The optimization can also be made faster using approximation. For example, optimization algorithm could be initialized using the results of spatial or temporal neighbors, and the parameters could be updated using a one-step Fisher’s scoring approximation. We would also like to explore more flexible methods for capturing the complicated wind motions, possibly a method where the window size is allowed to vary with location and time. This can ensure that our feature is always inside the frame of reference, resulting in more accurate estimates of the wind fields.


A. Standardization of data: We standardize the brightness temperature data as

where denotes the pixel-wise sample mean over time of at location , and denotes the corresponding smoothed standard deviation. The standard deviation is smoothed using the Gaussian kernel

where is the kernel bandwidth and denote the spatial locations over the entire region considered. The smoothing weights are defined as

That is, for standardizing the data we use

where is the sample standard deviation at location .

B. Smoothing the estimates: Let denote the estimated wind vector in using the proposed method. Also let the corresponding estimated variances be denoted by and where the suffixes ‘u’ and ‘v’ refer to the u- and v- component of the estimated wind vector. We smooth each component of the estimates using a weighted Gaussian smoothing filter, the weights being equal to the inverses of the corresponding variance estimates. Let denote the smoothed wind estimates where,


Supplementary Materials:

The supplementary material shows the 95% coverage probabilities for the u- and v- components of wind motion vector. (.pdf file)


The data analyzed in this paper were provided by Jessica Matthews of the North Carolina Institute for Climate Studies. This work was supported by the National Science Foundation (DMS - 1613219), the National Institutes of Health (R01ES027892) and King Abdullah University of Science and Technology (3800.2).


  • (1)
  • Ahmad et al. (2015) Ahmad, N., Suzuki, T. & Banno, M. (2015), ‘Analyses of shoreline retreat by peak storms using Hasaki Coast Japan data’, Procedia Engineering 116, 575–582.
  • Ailliot (2004) Ailliot, P. (2004), Modeles autorégressifsa changement de régimes Markovien-Applicationa la simulation du vent, PhD thesis, PhD Thesis. Université de Rennes 1.
  • Ailliot et al. (2006) Ailliot, P., Monbet, V. & Prevosto, M. (2006), ‘An autoregressive model with time-varying coefficients for wind fields’, Environmetrics 17(2), 107–117.
  • Bennett (1979) Bennett, R. J. (1979), ‘Spatial time series’, Pion, London 674.
  • Berkowitz & Steckelberg (2017) Berkowitz, B. & Steckelberg, A. (2017), ‘How Santa Ana winds spread wildfires’,
  • Boukhanovsky et al. (2003) Boukhanovsky, A. V., Krogstad, H. E., Lopatoukhin, L. J. & Rozhkov, V. A. (2003), Stochastic simulation of inhomogeneous metocean fields. part i: Annual variability, in ‘International Conference on Computational Science’, Springer, pp. 213–222.
  • Bras & Rodriguez-Iturbe (1985) Bras, R. L. & Rodriguez-Iturbe, I. (1985), Random functions and hydrology, Courier Corporation.
  • Brillinger (2001) Brillinger, D. R. (2001), Time series: Data analysis and theory, Vol. 36, Siam.
  • Brown et al. (1984) Brown, B. G., Katz, R. W. & Murphy, A. H. (1984), ‘Time series models to simulate and forecast wind speed and wind power’, Journal of Climate and Applied Meteorology 23(8), 1184–1195.
  • Castino et al. (1998) Castino, F., Festa, R. & Ratto, C. (1998), ‘Stochastic modelling of wind velocities time series’, Journal of Wind Engineering and Industrial Aerodynamics 74, 141–151.
  • Daniels et al. (2010) Daniels, J., Velden, C., Bresky, W., Wanzong, S. & Berger, H. (2010), GOES-R Advanced Baseline Imager (ABI) Algorithm Theoretical Basis Document for: Derived Motion Winds, University of Wisconsin–Madison.
  • De Luna & Genton (2005) De Luna, X. & Genton, M. G. (2005), ‘Predictive spatio-temporal models for spatially sparse enviromental data’, Statistica Sinica 15, 547–568.
  • ESA (2016) ESA, E. S. A. (2016), ‘Wind speed’,
  • Ester et al. (1996) Ester, M., Kriegel, H.-P., Sander, J., Xu, X. et al. (1996), A density-based algorithm for discovering clusters in large spatial databases with noise., in ‘Kdd’, Vol. 96, pp. 226–231.
  • Fuentes et al. (2008) Fuentes, M., Reich, B. J., Lee, G. et al. (2008), ‘Spatial–temporal mesoscale modeling of rainfall intensity using gage and radar data’, The Annals of Applied Statistics 2(4), 1148–1169.
  • Haas (1990) Haas, T. C. (1990), ‘Lognormal and moving window methods of estimating acid deposition’, Journal of the American Statistical Association 85(412), 950–963.
  • Haas (1995) Haas, T. C. (1995), ‘Local prediction of a spatio-temporal process with an application to wet sulfate deposition’, Journal of the American Statistical Association 90(432), 1189–1199.
  • Haslett & Raftery (1989) Haslett, J. & Raftery, A. E. (1989), ‘Space-time modelling with long-memory dependence: Assessing Ireland’s wind power resource’, Applied Statistics pp. 1–50.
  • Kim et al. (2014) Kim, T.-H., Yang, C.-S., Oh, J.-H. & Ouchi, K. (2014), ‘Analysis of the contribution of wind drift factor to oil slick movement under strong tidal condition: Hebei spirit oil spill case’, PLOS ONE 9(1), e87393.
  • Knapp & Wilkins (2018) Knapp, K. R. & Wilkins, S. L. (2018), ‘Gridded satellite (gridsat) GOES and CONUS data’, Earth System Science Data 10(3), 1417–1425.
  • Kyriakidis & Journel (1999) Kyriakidis, P. C. & Journel, A. G. (1999), ‘Geostatistical space–time models: A review’, Mathematical geology 31(6), 651–684.
  • Malmberg et al. (2005)

    Malmberg, A., Holst, U. & Holst, J. (2005), ‘Forecasting near-surface ocean winds with Kalman filter techniques’,

    Ocean Engineering 32(3-4), 273–291.
  • Modlin et al. (2012) Modlin, D., Fuentes, M. & Reich, B. J. (2012), ‘Circular conditional autoregressive modeling of vector fields’, Environmetrics 23(1), 46–53.
  • Monbet et al. (2007) Monbet, V., Ailliot, P. & Prevosto, M. (2007), ‘Survey of stochastic models for wind and sea state time series’, Probabilistic engineering mechanics 22(2), 113–126.
  • Priestley (1981) Priestley, M. B. (1981), Spectral analysis and time series, Academic Press London; New York.
  • Stein (2005) Stein, M. L. (2005), ‘Statistical methods for regular monitoring data’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67(5), 667–687.
  • Stein et al. (2013) Stein, M. L., Chen, J., Anitescu, M. et al. (2013), ‘Stochastic approximation of score functions for Gaussian processes’, The Annals of Applied Statistics 7(2), 1162–1191.
  • Tol (1997)

    Tol, R. S. (1997), ‘Autoregressive conditional heteroscedasticity in daily wind speed measurements’,

    Theoretical and Applied Climatology 56(1-2), 113–122.
  • Tomassini et al. (1999) Tomassini, M., Kelly, G. & Saunders, R. (1999), ‘Use and impact of satellite atmospheric motion winds on ecmwf analyses and forecasts’, Monthly Weather Review 127(6), 971–986.
  • WeatherOnline (2018) WeatherOnline (2018), ‘Calima’,