For example, the Santa Ana winds blow into California after scorching summers (Berkowitz & Steckelberg, 2017;Bopp & Shaby, 2017) and produce conditions that lead to large wildland fires. 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. There is a significant statistical literature on modeling winds from ground monitors (Brillinger, 2001;Haslett & Raftery, 1989;Priestley, 1981;Stein,2005b). For some applications, it is sufficient to analyze the evolution of wind at a fixed location, and several methods have been proposed for this scenario (Ailliot, 2004;Brown et al., 1984;Monbet et al., 2007;T o l ,1997). The wind at different locations can be utilized to model spatio-temporal dependencies (Bennett, 1979;Bras & Rodriguez-Iturbe, 1985;De Luna & Genton, 2005;Kyriakidis & Journel, 1999). 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 (2005b) 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. Sigrist et al. (2012) modeled precipitation with a nonseparable spatio-temporal model using an external wind vector. However, these wind observations have very sparse coverage. More recently, spatial modeling of winds were done at annual (Jeong et al., 2018), monthly (Jeong et al., 2019) and daily (Tagle et al., 2019;Tagle et al., 2020) temporal resolutions. Also, Chen et al. ( 2021) developed a non-Gaussian space-time model to assimilate observed data and simulations to predict future winds.

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, the satellites cannot monitor winds continuously in space and time, as they are in low earth orbit. Ground monitors can measure surface winds continuously in time but are sparse in space. 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.

On the other hand, geostationary weather satellites provide data from the surface and the atmosphere with a very high temporal resolution. The data comprise a series of images which essentially make them a "movie". While the satellites do not directly measure wind, they measure infrared channel radiations in terms of brightness temperature. These brightness temperature image sequences are used to infer wind estimates by tracking movements of atmospheric tracers such as clouds or moisture features over time. Figure 1 shows a sequence of brightness temperature images captured by the Geostationary Operational Environmental Satellite (GOES)-15, operated by the National Oceanic and Atmospheric Administration (NOAA). Wind data obtained from satellite images play a major role in data assimilation (Lahoz & Schneider, 2014). 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 spatially and temporally sparse. The DMW Algorithm also does not provide 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. Instead of modeling the covariance function globally, we borrow the idea of Nested Tracking from Daniels et al. (2010) by considering local data buffers and estimate wind vectors using maximum likelihood estimates under the assumption of local stationarity. Local likelihood estimation of covariance parameters has been studied extensively in Haas (1990); Haas (1995);Ver Hoef et al. (2004); Risser and Calder (2015); Kuusela and Stein (2018); Lenzi et al. (2020); Wiens et al. (2020); Castro-Camilo and Huser (2020). Estimating covariance parameters locally helps tackle the nonstationarity present in the data. One major advantage of our approach is that it circumvents the computational burden for large sample sizes. Moreover, the proposed algorithm can easily be run in parallel, which can lead to further computational improvements. Another major advantage of the approach over the DMW Algorithm is that it allows us to quantify uncertainties in the local estimates of the winds; this in turn facilitates our use of an inverse-variance-weighted spatial smoothing algorithm to borrow strength across space while down-weighting the most uncertain estimates. This also allows tradeoffs between bias and variance for reducing estimation error (Anderes & Stein, 2011). The estimated motion winds obtained from the proposed model should be assigned at the height of the features being tracked (subject to availability of altitude data), which will vary over space and time. See the Conclusions section for a detailed discussion of this issue.

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 5), 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 brightness temperature for 3 consecutive image times. The middle portion of the buffer is divided into smaller target scenes (or target windows), and each scene is analyzed to locate and select a set of suitable targets in the middle image. Daniels et al. (2010) also gives a description of the Nested Tracking Algorithm which involves nesting smaller target windows (usually of size 5 × 5) within a large target scene of size 15 × 15 pixels and getting every possible local motion vectors derived from each possible smaller window within a large target scene. The displacement vector between time points t and t + 1 is computed by minimizing the Sum of Squared Differences (SSD) criterion as

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

and is assigned as the DMW estimate at location x time point t in the buffer. Once every possible local motion vector 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 window 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 DMW Algorithm does not offer wind estimates at every space-time location as the data in the middle image has to satisfy a set of criteria to qualify as a suitable target scene (Daniels et al., 2010). Target scenes are processed to determine whether they (i) capture the intended target, (ii) contain sufficient contrast, and (iii) do not contain a mix of multiple potential targets. Target scenes that do not satisfy these criteria are deemed unsuitable for feature tracking and are rejected; no wind estimates are provided for those scenes. Also, quantifying uncertainties using the SSD criterion is hard because the subset of data used in the minimization criterion (1) depends on the parameter u, so likelihood ratio tests are not applicable. Finally, the highest resolution of the vector estimates generated for each target scene can be half-integers.

The proposed approach uses spatio-temporal covariance functions to track the wind. This requires us to consider asymmetric spatio-temporal covariance functions. Space-time asymmetries in covariance functions have been studied in Stein (2005a), Park and Fuentes (2006) and Huang and Sun (2019).

The space-time process Z(x, t) has asymmetric covariance if

Here (x, t) is a space-time location, with x =(x 1 , x 2 ) giving the spatial coordinates, and t giving the time. In most regions, winds flow in a consistent direction, and so changes in brightness temperature and other atmospheric variables at one location tend to precede similar changes down wind. For instance, if t 2 > t 1 and winds flow consistently from x to y,then we expect Cov{Z(x, t 1 ), Z(y, t 2 )} > Cov{Z(x, t 2 ), Z(y, t 1 )}.

We incorporate space-time asymmetries via a drift parameter. Suppose that Z 0 ∶ R 3 → R is a stationary, space-time symmetric process with covariance between Z 0 (x, t 1 ) and Z 0 (y, t 2 ) given by

where M 𝜈 is the isotropic Matérn (Stein, 2012) covariance function,

where 𝜅 𝜈 (⋅) is the modified Bessel function of the second kind of order 𝜈, Γ(⋅) is the gamma function, 𝜎 2 is the spatial process variance, and 𝜈 is the smoothness parameter which controls the mean square differentiability of the process. Also, a 11 and a 22 are the inverse spatial range parameters and a 33 is the inverse temporal range parameter. Z 0 is geometrically anisotropic in its spatial coordinates with anisotropy parameter a 12 , inverse spatial ranges a 11 and a 12 , and temporal range a 33 , but it is space-time symmetric because where the two-dimensional vector u =(u 1 , u 2 ) is the drift of the process over time and can be interpreted as the zonal and meridional components of wind. Intuitively, for a spatial process drifting from location x at time t 1 to location y at time t 2 , t 2 > t 1 , the optimal value of u can be thought of as the value that induces the highest correlation between the processes at (x, t 1 ) and (y, t 2 ). Then the covariance function of

which makes the covariance function in Equation (3) space-time asymmetric and stationary. This model not only incorporates wind motions via the drift parameter, but it also explicitly models how brightness temperature (or other suitable tracer field) naturally evolves over time, even in the absence of winds.

Assume the brightness temperature data Y (x, t) are normally distributed and have been standardized at each location (as described in the Appendix); denote the standardized data as Z(x, t). The mean and variance carry no information about the drift, and this step simplifies estimation of the correlation parameters. We do not specify a global covariance function for Z. Instead, we specify its local covariance with drift parameter u(x, t)∈R 2 and estimate the motion winds locally.

We define a target window as a square array of pixels

To estimate u(x, t), we assume that the process Z(x, t) is locally stationary in D(x, t) and that winds are smooth enough to be assumed constant in the scene, that is, u(x ′ , t ′ )≈u(x, t) for all (x ′ , t ′ )∈D(x, t). In other words, we approximate the local covariance function as

for (x 1 , t 1 ) and (x 2 , t 2 )∈D(x, t).

We use maximum likelihood estimation in local, moving windows to make inference about the spatially varying covariance parameters 𝜽(x, t)=(𝜎 2 (x, t), a 2 11 (x, t), a 2 22 (x, t), a 2 33 (x, t),𝜈(x, t), a 12 (x, t), u(x, t)). The local estimation approach deals with the nonstationarity present in the entire data. We use maximum likelihood estimation within D(x, t) to estimate 𝜽(x, t) ≡ 𝜽 D .IfZ D denotes the standardized data vector in the target window and Σ(𝜽 D ) denote the corresponding space-time covariance matrix with elements defined by (3), then the log-likelihood for

The estimated wind vectors are associated with the space-time location (x, t) at which the target window D was centered, denoted û(x, t). We also estimate the variances associated with the estimated wind vectors by computing the Fisher information at the maximum likelihood estimate, which gives a pointwise measure of uncertainty. We imitate the Nested Tracking approach and slide the target window across space and time, estimating wind vectors locally in space and time using the same optimization routine. In Section 4.1, we use the exact likelihood in (4) to estimate the wind vectors within a single target window. However, in Sections 4.2 and 5, we use a fast Vecchia approximation of the likelihood function (Guinness, 2018), implemented in the "GpGp" R package (Guinness & Katzfuss, 2018).

After obtaining the local estimates of u(x, t) for all (x, t), 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 window is an important tuning parameter. While implementing the method on real data sets, the window size is 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 the model parameters including window size on the performance of our method. We also analyze the efficacy of spatial smoothing of the local estimates. While implementing the methods in the real-data analysis in Section 5, the optimal window size is chosen using cross validation. This works well as wind fields are mostly smooth over a fairly small region and time frame.

In the first part of our simulation study, we determine the conditions under which the space-time drift model (STDM) performs well for estimating motion winds. For this purpose, we repeatedly simulate data sets within one particular target window (as opposed to scanning across a spatial domain) from a space-time Gaussian process with exponential covariance (obtained by substituting 𝜈 = 1∕2anda 12 = 0 in the expression for M 𝜈 in Equation ( 3)). We also implement a version of the DMW Algorithm and compare its performance with the STDM. To compare the two methods, accuracy for simulated data set i is measured by Vector Difference (Daniels et al., 2010) between the true (u 0 ) and estimated (û i ) wind vectors

and we report the mean and standard deviation of VD 1 , … , VD N over N = 100 data sets in Tables 123. First we consider repeated data simulations in a target window of size 7 × 7 × 3, generated independently from a Gaussian process with covariance function M 𝜈 with 𝜈 = 1∕2andu(x, t)=u 0 (constant true wind) using a fast Gaussian approximation algorithm (Guinness, 2018) which is implemented using the "GpGp" package (Guinness & Katzfuss, 2018) in R. We assume 1∕a 11 = 1∕a 22 = 𝛼 1 and the true spatial range parameter, 𝛼 1 is chosen to be either 1, 2, 4, or 8. The true temporal range parameter 1∕a 33 = 𝛼 2 is chosen to be either 1, 2, 3, or 4. We also take two different values of the true wind vector, namely u 0 =(1, 2) T and (3, 5) T , which signify respectively slow and fast wind vectors. Except 𝜈, which is fixed at 𝜈 = 1∕2, all covariance parameters are updated simultaneously during optimization. Table 1 shows the performance of the two methods for the two wind vectors based on N = 100 simulations.

STDM does a better job in estimating moderately small and large wind vectors for the 7 × 7 window size compared to DMW Algorithm. 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. In particular, windows of size 7 × 7 are not adequate to contain the features being tracked in frame for a wind vector of (3, 5) T which is reflected in the high MVD and SD values for STDM (see Table 1, top panel).

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 Algorithm are at most half integers and they are limited to the size of the larger search window. That is, while estimating motion vectors, the smaller central target scene (3 × 3) 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 also be attributed to the window size relative to the magnitude of the wind vector.

To examine the effect of window size, we now perform the simulations again with window sizes of 11 × 11 and 15 × 15 respectively, under the same covariance parameter settings as described earlier. The panels of Tables 2 and 3 portray a clear picture of the effect of window size on the estimation of winds using STDM. The estimation of both wind vectors improve as we increase window size from 7 × 7to11× 11 (see Table 2) and then to 15 × 15 (see Table 3). for true wind vectors u 0 =(1, 2) T (left panel) and u 0 =(3, 5) T (right panel) based on data window of size 7 × 7; 𝛼 1 and 𝛼 2 denote the true spatial and temporal range respectively. MVD and SD are measured in the number of pixels. In this part of the simulation study, the true wind has been chosen in a single window. This is analogous to constant wind vectors over a spatial domain. The next section of the simulation studies consider nonconstant wind fields and investigate the performance of the proposed method (including spatial smoothing) for different window sizes when data are generated from nonstationary/non-Gaussian processes. An additional simulation setting can be found in the Appendix.

We now conduct another simulation study which examines the performance of the competing methods and looks into the efficacy of spatial smoothing of wind estimates under three different data generation schemes. Unlike the simulation setting in Section 4.1, we now generate data sets on a 100 × 100 unit square at three consecutive time points. While generating the data, we consider rotational winds, which are obtained by rotating each location counter-clockwise by 6 • around the origin at each time point. If x denotes a spatial location at time t and x ′ denotes the same at time (t + 1),then, for t = 1, 2and𝜃 = 6𝜋∕180. Then

These winds are nonuniform over space and emulate big tropical storms which rotate counter-clockwise in the Northern Hemisphere and clockwise in the Southern Hemisphere due to the Coriolis effect. This simulation setting reflects realistic scenarios and compensate for the lack of availability of data sets with ground truth for wind measurements.

Under the first data generation scheme, the shifted locations at three consecutive time points are used to generate data sets from an approximate Gaussian process using the "GpGp" package (Guinness & Katzfuss, 2018) in R. The covariance function used for generating the data is the exponential space-time covariance (M 𝜈 with 𝜈 = 1∕2) with true spatial range 𝛼 1 = 0.075 (4 pixels), true temporal range 𝛼 2 = 5 and a nugget variance of 0.1. The second data generation model is a Gaussian process model with nonstationary Matérn covariance function (Paciorek & Schervish, 2006) with nugget. In particular, we now assume that the spatial range parameter in the covariance function is spatially varying. The nonstationary covariance function has been defined at the shifted locations, so as to incorporate the rotational shift in the generated data. Analytically, the space-time version of the nonstationary Matérn covariance function can be written as

where Σ=Σ(x)=

x 1 5 I 2 and Σ ′ =Σ(x ′ ) are kernel covariances which evolve over the spatial domain to produce nonstationary covariances. Here, I 2 denotes the identity matrix. For our simulation, we choose 𝜎 2 = 2,𝜈 = 1.5a n d 𝜏 2 = 0.1. Finally, for our third simulation setting, we generate data from a non-Gaussian process, in particular, a Gamma process. For this purpose, we first generate Gaussian data Z(x ′ , t) at the shifted locations x ′ for three consecutive time points, and then set

where F is the CDF of a Gamma distribution with parameters (a, b) and Φ is the standard normal CDF. We choose a = b = 2 in our simulation. Since we would never know if the data were exactly Gaussian, in practice, we often fit a Gaussian process to non-Gaussian data. Thus, the goal of this simulation setting is to evaluate the efficacy of the proposed local model under slight deviations from normality.

From each of the three data generating schemes, we simulate N = 50 data sets and fit the proposed space-time drift model locally. For data generated under the exponential space-time covariance model, we estimate all covariance parameters except the smoothness parameter, which is fixed at 𝜈 = 0.5. For the other two settings, all covariance parameters are simultaneously estimated. We estimate the parameters at 75 × 75 locations at the center of the grid. This reduces boundary effects during local estimation and allows us to compare the performance of our method for multiple window sizes. The wind estimates are spatially smoothed using a weighted Gaussian kernel. We also estimate the winds using DMW Algorithm under similar settings, and the DMW estimates are also smoothed for fair comparison. For both methods, the optimal smoothing bandwidth has been chosen using cross-validation. Table 4 compares the average performance of the algorithms and their smoothed versions over N = 50 data sets, simulated from the three different settings. From the results, it can be seen that the STDM outperforms the DMW Algorithm for all window sizes, expect for a large window size like 25 × 25 under the non-Gaussian setting for which the performances are comparable. Also, under every setting, smoothing the estimates gives us a much better representation of the wind fields compared to the raw estimates and STDM outperforms DMW Algorithm once smoothing is implemented. Also, Figures 234show the respective true wind fields (in red) along with the raw and smoothed estimates (in blue) from both methods for a single representative data set from each of the three settings respectively, obtained using windows of size 25 × 25. The wind vectors are plotted for a subset of the 75 2 locations in the spatial domain for clear visualization of the true and estimated fields. From the figures, it can be seen that the space-time drift model, fitted locally, efficiently captures the overall nonstationarity in the spatial fields and performs well under deviations from normality. The size of the target window plays a crucial role while estimating motion winds. For example, while estimating the rotational wind field using 15 × 15 (or smaller) windows, some of the target features move out of scene because of the nonuniformity in winds across space. In that scenario, irrespective of the data generating scenario, both methods provide poor estimates of the wind field. Using larger target windows, such as windows of size 25 × 25, fixes the problem to some extent. This corresponds to results from experimental studies carried out by special GOES-10 rapid scan imagery which showed that a higher number of high quality winds could be derived with decreasing time intervals and increasing spatial resolution (Daniels et al., 2010). For rotational wind field, larger window size improves estimation of the wind fields. This is because the average of the rotational wind field over a window is equal to the wind vector at the center of the window. Thus, a major takeaway from the simulation studies is that while estimating nonuniform wind fields, selecting the window size requires balancing a trade-off 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. Also, there is computational trade-off, since larger window sizes will lead to heavier computations. These motivate us to choose the window size using cross validation in the real-data analysis.

APPLICATION TO GOES-15 DATA

The Geostationary Operational Environmental Satellite (GOES), operated by National Oceanic and Atmospheric Administration (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 data set 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 (Kleiber et al., 2013;Paciorek & Schervish, 2006). Central Colorado is home to the Rocky Mountains, which reaches elevations of over 4000 m. To the east of the Rockies, the terrain gradually transitions into the Great Plains, dipping to a minimum elevation of approximately 1000 m. Thus, meteorological processes evolving across the domain are likely nonstationary and cannot be modeled using simple isotropic covariance functions. Moreover, several studies have concluded that winds travel faster, and often turn into windstorms, as they flow along the downward slope of the Continental Divide (Mercer et al., Also, Colorado has historically experienced windstorms during early January, such as the large-magnitude wind gusts in 2007 (Mercer et al., 2008). Therefore, the chosen spatio-temporal domain ensures that the spatial process modeled is nonstationary and the wind fields to be estimated are nonuniform, which helps best illustrate the proposed method.

Figure 5 shows the data in two representative data buffers centered at 00:15 am on January 03, 2015 (t = 193; top panel) and at 08:45 am on January 04, 2015 (t = 322; bottom panel) respectively. These two time points represent wind movements across the state at night and in the morning respectively. Lower values of brightness temperature indicate cloud cover, whereas high brightness temperature values suggest clear skies over the region.

The first step involves standardizing the brightness temperature data using the pixel-wise sample mean and standard deviation over time. This removes the constant background image and the effect of low cloud cover over the region which can affect the local estimation of wind vectors. To smooth the standard deviation map, we use a Gaussian kernel with smoothing parameter 𝜆 = 3 pixels (i.e., 12 km). We fit a Gaussian process with the full covariance function in (3) to the standardized data within each target window using the 'matern_anisotropic3D_alt' covariance function in the "GpGp" package in R. As mentioned earlier, the covariance function in (3) includes an additional parameter 12 to allow for anisotropy in the spatial covariances. All covariance parameters, including the smoothness parameter, are estimated locally within each window. We estimate the parameters at locations in the center of our spatial domain, so that the target windows are always a square around the target locations. We use cross validation to determine the appropriate size of target windows. Since no direct measurements of winds are available, we compare window sizes indirectly based on 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 all subsequent time points. Table 5 gives the Mean Squared Prediction Error (MSPE) and the Continuous Rank Probability Score (CRPS) (discussed in Section 5.3) for different window sizes and the corresponding average computation time (in seconds) to estimate the wind vector at one pixel using 3 consecutive time steps. All computations are done on a PC with Intel Core i7-9750H CPU at 2.60GHz with 32 GB of RAM. Based on a lower cross validation MSPE, comparable CRPS and a much lower computation time, we select sub regions of size 25 × 25 pixels (i.e., 100 km × 100 km) and estimate the wind vectors, along with all other covariance parameters, locally at each spatial location and at each time using maximum likelihood estimates. We also estimate the variances associated with the estimates by computing the Fisher information at the MLE. Figure 6 shows the uncertainty associated with the estimated wind components as given by the estimated standard deviations in the log scale, for the two data buffers mentioned above.

The estimated variances are high at some locations, especially where it is difficult to distinguish clouds from the background terrain and track their movements across the region. 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.

Similar to the DMW Algorithm, we also impose some simple criteria to check whether a target window D(x, t) is suitable for tracking and estimating motion winds. We look at the mean of the residual variances from regressing Z D (⋅, t) on Z D (⋅, t -1) and Z D (⋅, t + 1) on Z D (⋅, t).W es a yt h a tD(x, t) is discarded if the mean residual variance is low and the brightness temperature at the location, Y (x, t) is high. This signifies no cloud cover or negligible movement across time within the target window. For our data, we reject D(x, t) if the mean residual variance is in the lower 10 percentile of all mean residual variances and the brightness temperature Y (x, t) is in its upper 10 percentile, which is, Y (x, t) > 260 Kelvin. Figure 7 shows the raw and smoothed wind estimates obtained from the proposed STDM. The wind vectors have been plotted at a subset of spatial locations for better visualization.

For 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 DMW Algorithm are plotted at a subset of locations in Figure 8. The results once again highlight the limitations of the DMW Algorithm output. While the space-time drift model is able to capture the consistent winds over Colorado, the DMW Algorithm does not identify the pattern. As a result, the wind fields from the DMW Algorithm look jumbled after smoothing.

For this dataset, reference wind fields are not available. Hence, we compare the two methods based on Mean Squared Prediction Error (MSPE) and Continuous Rank Probability Score (CRPS; see (Gneiting & Katzfuss, 2014; Gneiting & Raftery, 2007) for details) while predicting standardized brightness temperature fields. To predict Z at a spatial location x time t, we consider Z D (⋅, t -1), the standardized data in D(x, t -1) and estimated winds at time t -2. This is because the estimated winds at time t -2 uses data from times t -3, t -2andt -1. The predicted standardized temperature is calculated as the mean of the Gaussian conditional distribution of Z(s, t) given Z D (⋅, t -1) under the stationary Matérn space-time drift model with estimated drift parameter û(s, t -2). We predict standardized brightness temperature using the raw and smoothed versions of the wind estimates from both algorithms. We also consider a naive approach of predicting the brightness temperature fields, where the data at time -1 is considered to be the predicted standardized brightness temperature fields at time t. We call it the persistence prediction and the performances of the two methods are assessed relative to the baseline MSPE. To calculate the CRPS, we assume the associated predictive distribution of standardized brightness temperature is well approximated by a Gaussian distribution with mean centered at the predicted value. We are also able to compute CRPS for predictions based on motion winds derived from DMW Algorithm by predicting standardized brightness temperature using the proposed probabilistic model. While making these predictions for Algorithm, we use all estimated covariance parameters from the STDM, except the drift parameters. Table 6 compares the raw and smoothed wind estimates in terms of MSPE and CRPS for 4 consecutive time points corresponding to the two data buffers. From Table 6, we also see a moderate gain in prediction accuracy for these time points while predicting using the smoothed wind estimates, suggesting that local smoothing is preferred. Figure 9 brightness temperature at subset of spatial locations and at t = 196 and t = 325 respectively using the smoothed estimated winds from STDM. It can be seen that our model captures the main features in the brightness temperature fields over Colorado. Our model also captures the wind movement over the region since we can see similar patterns of features across the region as compared to the ones tracked along in the original images (see Figure

Wind is one of the most important atmospheric variables; it has large impact on local weather and hence, studying winds is essential. In this paper, we propose a spatio-temporal model to derive atmospheric motion winds data from a sequence of high-resolution of images over time collected by geostationary satellites. In particular, we use a sequence of brightness temperature images taken by the GOES series of the NOAA meteorological satellites and estimate atmospheric wind speed and direction. We model brightness temperature using a space-time drift process and propose 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 paradigm of the Nested Tracking implementation of the Derived Motion Winds Algorithm, we propose local likelihood estimation of the covariance parameters. 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. We also compare the performance of our method with the DMW Algorithm. Section 4.1 details extensive simulation studies that outlines conditions under which our model performs well. Based 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 in Section 4.2 show that our method performs well even when we have nonstationary and non-Gaussian spatial fields. The simulations also show that spatial smoothing provides better estimates of the wind fields. The simulations highlight a major drawback of the DMW Algorithm. Due to the design of the DMW Algorithm, the local DMW Algorithm 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 limitation of the DMW Algorithm stems from the fact that the algorithm requires the original and shifted brightness temperature data to be at the same locations, in order to calculate the sum of squares criterion between the two images. However, this is not a requirement for our method, since we estimate the drift by maximizing a multivariate normal likelihood over a model that involves a continuous parameter u, rather than relying on the sum of squared differences between the original and shifted images. 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 from a small window. apply our method on brightness temperature data over 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, showing that our method is capable of quantifying uncertainties associated with the estimation. Since the uncertainties in estimated wind speed are correlated over space and time and hence hard to estimate, we provide point-wise uncertainties. We also provide smoothed maps of estimated wind fields over Colorado. We use cross-validation to choose the smoothing bandwidth thereby implementing a data-driven method for how much smoothing to do. Since zonal (horizontal) and meridional (vertical) wind components are usually studied separately in various applications, we smooth the two components of the wind vector separately. However, one can smooth both wind components together using a bivariate Gaussian smoothing filter, thereby using the correlation between the estimated wind components, along with the estimated variances. We predict brightness temperature fields using our model and compare the raw and smoothed estimates of wind fields based on Mean Squared Prediction Error (MSPE) and Continuous Rank Probability Score (CRPS). We also compare the performance of our proposed method and the DMW Algorithm. We have shown that our method outperforms the DMW Algorithm with respect to MSPE and CRPS. We argue that smoothing the estimated wind fields give more reliable estimates. We also see that we capture the main features in the brightness temperature maps through our prediction, including the drift across the region.

The estimated motion winds obtained from the DMW Algorithm are used as one of the input variables in the Multi-platform Tropical Cyclone Surface Wind Analysis (MTCSWA) model (Knaff et al., 2011), whose main purpose is to create an algorithm package to operationally generate an estimation of surface wind fields around active tropical cyclones. It has also been seen that, considering uncertainties of parameters usually results in improved projections from these large-scale scientific models (Sharma et al., 2023). As such, more accurate estimates of motion winds along with their estimated uncertainties obtained from the proposed method, can potentially improve the forecasts obtained from the MTCSWA model.

Wind speeds tend to increase with height above ground and as a result, wind speeds are usually measured at varying altitudes. According to Daniels et al. (2010), the GOES-R ABI cloud height algorithm is usually used to assign a representative height to each cloudy target. This algorithm uses National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) forecast temperature profiles and cloud-top pressure as ancillary data to assign a representative height to clouds in a target scene. However, since neither the output from the GOES-R ABI cloud height algorithm nor the ancillary data sets were available to us, we were unable to assign appropriate altitude levels to the estimated motion winds. If the heights were available, though, we would simply assign our estimated winds to the inferred top height of the clouds. Integrating these ancillary data sets to estimate the altitudes of the motion winds is nontrivial, and can be explored in a future project. The unavailability of data on cloud heights does make comparisons of estimated winds against like wind observations more difficult. Hence, we devised a validation scheme based on predicting future brightness temperatures based on the estimated wind fields.

The DMW Algorithm rejects target scenes that do not satisfy certain criteria and does not provide any wind estimates for those scenes. In our data application, we have tried to emulate the general setting of the Nested Tracking implementation of the DMW Algorithm, and have applied one of their target selection criterion on our data set to make the algorithms as comparable as possible. If the goal of the study is simply to report estimated winds when suitable tracer fields are available, then we think that the mechanism for missingness is not important. However, if one wants to use the estimated winds to infer long-term wind climatologies, then one should investigate the missingness mechanism. Both our method and the competitor are more likely to be missing on days without suitable tracers, and so if there is an association between cloud cover and wind then this is a source of informative missingness that should be taken into account inputting missing values. We expect such relationships to be weak, but further study is warranted for long-term climatological inference.

One major challenge 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 parallel 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, the 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. in Figure A1.FromFigureA1, it can be seen that the estimates of the smaller wind vectors on both sides of the central meridian become less accurate with increasing window size. This accounts for the increase of the MVD and the associated SD for window sizes larger than 25 × 25.