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 oilspills (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 largescale 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 & RodriguezIturbe, 1985; Kyriakidis & Journel, 1999; De Luna & Genton, 2005). Also, Boukhanovsky et al. (2003), Malmberg et al. (2005) and Ailliot et al. (2006) proposed autoregressive spacetime models to describe the evolution of winds. Stein (2005) proposed a spectralintime modeling approach to describe the spacetime 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 MediumRange 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 spacetime 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 2dimensional 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 aswhere 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 twodimensional 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 densitybased 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 halfintegers.
3 Modelbased wind estimation
3.1 Spacetime drift models
The proposed approach uses spatiotemporal covariance functions to track the wind. This requires us to consider asymmetric spatiotemporal covariance functions. The spacetime process has asymmetric covariance if
(1) 
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 spacetime asymmetries via a drift parameter. Suppose that is a stationary, spacetime symmetric process with covariance function
and let
Then the covariance function of is
(2)  
which is spacetime 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
(3) 
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 spacetime covariance matrix with elements defined by (3), then the loglikelihood for given is
(4) 
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 realdata 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 spacetime 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
1.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 
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 vcomponent 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 vcomponents 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 
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 
5 Application to GOES15 data
5.1 GOES15 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 GOES8 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 GOES15 satellite. Launched in March 2010, GOES15 is positioned at the GOESWest 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 pixelwise 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 GOES15 data can be obtained from NOAA OneStop at https://data.noaa.gov/onestop/#/collections.
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 southwest to northeast. 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.
5.2 Estimation using GOES15 data
The first step involves standardizing the brightness temperature data using the pixelwise 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 
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.
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 spacetime 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.
Method  

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 
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.
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 highresolution 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 GOESR 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 halfintegers 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 GOES15 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 onestep 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.
APPENDIX
A. Standardization of data: We standardize the brightness temperature data as
where denotes the pixelwise 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 MATERIAL
 Supplementary Materials:

The supplementary material shows the 95% coverage probabilities for the u and v components of wind motion vector. (.pdf file)
ACKNOWLEDGEMENTS
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).
References
 (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 MarkovienApplicationa 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 timevarying 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’, https://www.washingtonpost.com/graphics/2017/national/santaanafires/.
 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 & RodriguezIturbe (1985) Bras, R. L. & RodriguezIturbe, 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), GOESR 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 spatiotemporal models for spatially sparse enviromental data’, Statistica Sinica 15, 547–568.
 ESA (2016) ESA, E. S. A. (2016), ‘Wind speed’, https://sentinel.esa.int/web/sentinel/userguides/sentinel3altimetry/overview/geophysicalmeasurements/windspeed.
 Ester et al. (1996) Ester, M., Kriegel, H.P., Sander, J., Xu, X. et al. (1996), A densitybased 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 spatiotemporal 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), ‘Spacetime modelling with longmemory 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 nearsurface ocean winds with Kalman filter techniques’,
Ocean Engineering 32(34), 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(12), 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’, https://www.weatheronline.co.uk/reports/wind/TheCalima.htm.
Comments
There are no comments yet.