This post is an account of my recent experience with Google Earth Pro (GPE), a free online tool meant to make visualization of geographical information intuitive and accessible. This account contains tidbits that I thought others might find useful, but it is not meant to be an unbiased or comprehensive resource.
My goal in using GPE was to set up a short video tour of the reservoirs on the main stem of the Upper Snake River, which has its upper reaches around the Teton Range (WY, USA) before going through much of southern Idaho. Here is the video, produced with instructions from this post from the blog:
The cool
GPE is indeed an intuitive way for people that have only a limited experience of geographic information systems (GIS) to put together nice-looking videos in a limited amount of time. It relies on the increasingly large amounts of geographic data available online. For instance, polygons made of points with geographical coordinates, used to delineate political boundaries or catchments among others, can be found increasingly easily, be it under formats specific to GPE (KML, KMZ), or traditional GIS format such as SHP (shapefile format, which can be imported to GPE with… “File => Import”, I told you it was intuitive). This capability to find the right data could be improved by the recent launch of a dataset search engine by Google.
This video, after several tests, superposed three layers on the satellite images of the land surface:
- Pins of the dams’ locations. Please refer to this post to learn all you need to know about setting and customizing such new placemarks.
- The network of all streams with an average flow of more than 10 m3/s in the Columbia River basin. This is mainly for easy visualization of where the Snake River and its affluents are. That data was obtained from this useful USGS resource, with all the major basins in the US having data in GPE-ready KML format (basin boundaries, tributaries with their average flows, landuse, etc.).
- An invisible layer contains the shapes of the reservoirs’ lakes, in order to zoom in on the reservoir and not the dam itself. I got the shapes of waterbodies in the Upper Snake River basin from an amazing resource for water resources practitioners (see “The limits” below to see how to access this resource, but also to understand why it must be handled with care). Hint: since the list of waterbodies is dauntingly large and includes every little puddle in the area, I had to zoom in to the desired reservoir so as to only select its shape, by doing “File => Import”, then when prompted, by choosing to only import the features in my line of sight.
The trade-offs
The trade-offs when doing this kind of short intro video to your study area, are between video quality and required memory. This is especially true if the video is of a rather large area, as the GPE satellite images embed a level of detail that is amazing at the scale of that basin. Therefore, when creating my video, I had to resort to several tries before getting a correct quality, and this video eats up 500Mo of disk space for 43 seconds (!!!!). Any video with lower resolution just looked downright disgusting, whereas any larger video may have serious problem running on your laptop. Think of it as a Goldilocks zone that you have to find if you want to have a video you can embed in an oral presentation.
(NB: to embed this video in a presentation in a fullproof kind of way, the easiest is to embed a link from the PPT file to the video in the same folder).
The limits
The limits of Google Earth are with its plain refusal to display features it deems too big. To understand this, let us look at this USGS webpage where a sizable of hydrological information can be downloaded. In particular it is possible to download data exactly for the upper Snake area (Subregion 1704 on the picture below). This data can include waterbodies as discussed above (in “the cool”) and that is contained in the NHD data ticket in the picture. By ticking “Watershed Boundary Dataset (WBD)”, one can also download basin and subbasin boundaries.
Why then did I not represent Upper Snake boundaries in my video? Well, the Upper Snake polygon, under SHP format at least, is too big so GPE just… refused to display it. I tried to represent the HU-6 and HU-8 subbasins (smaller than the HU-4-sized Upper Snake under this classification system, and:
- The HU-6 subbasin just divides the basin into “headwaters” (smaller part upstream of Palisades) and “the rest”; the former feature was smaller and GPE plotted it, but it did not plot the latter.
- The HU-8 subbasins all are individually smaller features that GPE accepts to plot. So I could have plotted the watershed boundaries by plotting ALL of the HU-8 subbasins, but spoiler alert: this looked horrendous.
Takeaway: GPE only displays features that individually have limited size. So maybe using a more powerful and recent desktop computer than the one I had would have done the trick… but be aware that there will always be a limit. Also note that assuming that GPE is just taking its time loading the data, and that staring blankly at the screen in the meantime is not going to be too painful, is not a good idea. Take a nap instead, or do something else, then admit that the damn thing just did not display.
A solution for this, of course, is to load the larger features on GIS software and making them lighter by eliminating points in the polygon without altering the shape… but that kinda beats the purpose of GPE, which is to avoid having to become a GIS geek, doesn’t it?
for which we are trying to estimate 
with Bayesian OLR is
, 
being the independent variables and 
the dependent variable, 
is a vector of model parameters and 
is the model error variance. The unusual notation for the normal distribution 
means “a normal distribution of 
and variance 
,” where 
is the number of points in the data set. The estimated variance of 
, 
can be calculated as
and 
are the mean and the covariance of the 
with an error variance of 
), replace 
and 
, e.g. 
, to add features to 
, we would instead have that
and 
. One problem with this approach is that as we add more terms (also called features) to function 
to better capture non-linearities, the the dimension of the matrix we will have to invert increases. It is also not always obvious what features we should add to our data. Both problems can be handled with Gaussian Processes.
we have to invert (remember the dimensions of such matrix equals the number of features in 
— 
is just another point 
either in the data set or for which we want to estimate 
will be called a covariance or kernel function. Since the prior’s covariance 
. Using 
and assuming a prior of 
, 
, 
and 
. Now the features of the data are absorbed within the inner-products between observed 
‘s and or for 
, is that 
does not have to represent simply a matrix-matrix or vector-matrix multiplication of 
, which will necessarily return a positive semi-definite matrix and therefore be a valid covariance matrix. There are several kernels (or kernel functions or kernel shapes) available, each accounting in a different way for the nearness or similarity (at least in Gaussian Processes) between values of 
:
. This kernel is used when trying to fit a line going through the origin.
, where d is the dimension of the polynomial function. This kernel is used when trying to fit a polynomial function (such as the fourth order polynomial in the blog post about Bayesian OLR), and
, where 
is a kernel parameter denoting the characteristic length scale of the function to be approximated. Using the RBF kernel is equivalent to adding an infinite (!) number of features the data, meaning 
.
at discrete points, with a collection of 
for all 
. After all, don’t we want a function more often than not solely for the purpose of getting point estimates and associated errors for given values of 
. The mean of all possible functions at 
in 
will be for each pair 
based on the value of a distance metric between 
tend to be similar and vice-versa. We can therefore turn the sampled functions 
into whatever functional form we want by choosing the appropriate covariance matrix (kernel) — e.g. a linear kernel will return a linear model, a polynomial kernel will return a polynomial function, and an RBF kernel will return a function with an all linear and an infinite number of non-linear terms. A Gaussian Process can then be written as
is the mean function (e.g. a horizontal line along the x-axis at y = 0 would mean that ![\mu(\boldsymbol{x}) = [0,..,0]^T](https://s0.wp.com/latex.php?latex=%5Cmu%28%5Cboldsymbol%7Bx%7D%29+%3D+%5B0%2C..%2C0%5D%5ET&bg=ffffff&fg=444444&s=0&c=20201002)
) and 
along the x-axis from x=-5 to x=5 (
) using a Gaussian Process. Assume the parameters of our Gaussian Process from which we will sample our functions are mean 0 and that it uses an RBF kernel for its covariance matrix with parameter 
(not be confused with standard deviation), meaning that ![k_{i,i}=1.0 ,\:\forall i \in [0, 199]](https://s0.wp.com/latex.php?latex=k_%7Bi%2Ci%7D%3D1.0+%2C%5C%3A%5Cforall+i+%5Cin+%5B0%2C+199%5D&bg=ffffff&fg=444444&s=0&c=20201002)
and ![k_{i,j}=e^{-\frac{||x_i-x_j||^2}{2\cdot 2^2}},\: \forall i,j \in [0, 199]](https://s0.wp.com/latex.php?latex=k_%7Bi%2Cj%7D%3De%5E%7B-%5Cfrac%7B%7C%7Cx_i-x_j%7C%7C%5E2%7D%7B2%5Ccdot+2%5E2%7D%7D%2C%5C%3A+%5Cforall+i%2Cj+%5Cin+%5B0%2C+199%5D&bg=ffffff&fg=444444&s=0&c=20201002)
, or![\begin{aligned} p(\boldsymbol{f}_*|\boldsymbol{X_*})&=\mathcal{N}(\boldsymbol{f}|\boldsymbol{\mu},\boldsymbol{K_{**}})\\ &=\mathcal{N}\left(\boldsymbol{f}_*\Bigg|[\mu_{*0},..,\mu_{*199}], \begin{bmatrix} k_{*0,0} & k_{*0,1} & \cdots & k_{*0,199} \\ k_{*1,0} & k_{*1, 1} & \cdots & k_{*1,199}\\ \vdots & \vdots & \ddots & \vdots \\ k_{*199,0} & k_{*199,1} & \cdots & k_{*199,199} \end{bmatrix} \right)\\ &=\mathcal{N}\left(\boldsymbol{f}_*\Bigg|[0,..,0], \begin{bmatrix} 1 & 0.9997 & \cdots & 4.2\cdot 10^{-6} \\ 0.9997 & 1 & \cdots & 4.8\cdot 10^{-6}\\ \vdots & \vdots & \ddots & \vdots \\ 4.2\cdot 10^{-6} & 4.8\cdot 10^{-6} & \cdots & 1 \end{bmatrix} \right) \end{aligned}](https://s0.wp.com/latex.php?latex=%5Cbegin%7Baligned%7D+p%28%5Cboldsymbol%7Bf%7D_%2A%7C%5Cboldsymbol%7BX_%2A%7D%29%26%3D%5Cmathcal%7BN%7D%28%5Cboldsymbol%7Bf%7D%7C%5Cboldsymbol%7B%5Cmu%7D%2C%5Cboldsymbol%7BK_%7B%2A%2A%7D%7D%29%5C%5C+%26%3D%5Cmathcal%7BN%7D%5Cleft%28%5Cboldsymbol%7Bf%7D_%2A%5CBigg%7C%5B%5Cmu_%7B%2A0%7D%2C..%2C%5Cmu_%7B%2A199%7D%5D%2C+%5Cbegin%7Bbmatrix%7D+k_%7B%2A0%2C0%7D+%26+k_%7B%2A0%2C1%7D+%26+%5Ccdots+%26+k_%7B%2A0%2C199%7D+%5C%5C+k_%7B%2A1%2C0%7D+%26+k_%7B%2A1%2C+1%7D+%26+%5Ccdots+%26+k_%7B%2A1%2C199%7D%5C%5C+%5Cvdots+%26+%5Cvdots+%26+%5Cddots+%26+%5Cvdots+%5C%5C+k_%7B%2A199%2C0%7D+%26+k_%7B%2A199%2C1%7D+%26+%5Ccdots+%26+k_%7B%2A199%2C199%7D+%5Cend%7Bbmatrix%7D+%5Cright%29%5C%5C+%26%3D%5Cmathcal%7BN%7D%5Cleft%28%5Cboldsymbol%7Bf%7D_%2A%5CBigg%7C%5B0%2C..%2C0%5D%2C+%5Cbegin%7Bbmatrix%7D+1+%26+0.9997+%26+%5Ccdots+%26+4.2%5Ccdot+10%5E%7B-6%7D+%5C%5C+0.9997+%26+1+%26+%5Ccdots+%26+4.8%5Ccdot+10%5E%7B-6%7D%5C%5C+%5Cvdots+%26+%5Cvdots+%26+%5Cddots+%26+%5Cvdots+%5C%5C+4.2%5Ccdot+10%5E%7B-6%7D%C2%A0%26%C2%A04.8%5Ccdot+10%5E%7B-6%7D+%26+%5Ccdots+%26+1+%5Cend%7Bbmatrix%7D+%5Cright%29+%5Cend%7Baligned%7D&bg=ffffff&fg=444444&s=0&c=20201002)
![\boldsymbol{f}_* = [y_0, ..., y_{199}]](https://s0.wp.com/latex.php?latex=%5Cboldsymbol%7Bf%7D_%2A+%3D+%5By_0%2C+...%2C+y_%7B199%7D%5D&bg=ffffff&fg=444444&s=0&c=20201002)
. A draw of 3 such sample functions would look generally like the following:
. All we have to do is to add random variables corresponding to our data and condition the resulting distribution on them, which means accounting only for samples of functions that exactly go through (or given) our data points. The resulting distribution would then look like 
. After adding our data 
, or 
, our distribution would look like
from which we want to get point estimates and corresponding error standard deviations. The covariance matrix above accounts for correlations of all observations with all other observations, and correlations of all observations with all points to be predicted. We now just have to 
, 
, 
, and 
is a vector containing the elements in the diagonal of the covariance matrix 
— the variances 
of each random variable ![\boldsymbol{\mu} = \boldsymbol{\mu}_* = [0, ... , 0]^T](https://s0.wp.com/latex.php?latex=%5Cboldsymbol%7B%5Cmu%7D+%3D+%5Cboldsymbol%7B%5Cmu%7D_%2A+%3D+%5B0%2C+...+%2C+0%5D%5ET&bg=ffffff&fg=444444&s=0&c=20201002)
, the mean of the Gaussian Process would be 
. A plot of 
estimated for 200 values of 
), regressed on 5 random data points 
should look generally like
of the RBF kernel, which determines the correlation length scale of the fitted model and therefore how fast with increasing separation distance between data points the model will default to the prior (here, with 
). The figure below exemplifies this effect
, we’d do something like:
and variance 
for each of the regression model parameters 
is defined as
is the mean of the model error (normally 
), 
. Digression: in Generalized Linear Regression (GLR), observations may not be independent and may have different variances (heteroscedasticity), so the covariance matrix may have off-diagonal terms and the diagonal terms may not be equal.
for this function, but do not confuse it with a normal distribution. This procedure is called basis function expansion and has a likelihood function of the form
, called the parameter posterior distribution, is calculated by multiplying the likelihood function used in MLE by a prior distribution for the parameters. A prior distribution is assumed from knowledge prior to analyzing new data. In Bayesian Linear Regression, a Gaussian distribution is commonly assumed for the parameters. For example, the prior on 
for algebraic simplicity. The parameter posterior distribution assuming a known 
, given 
is a constant depending only on the data.
for a data set 
for data set 
, with 
denoting a new point that is not in 
![\phi(x) = [1, x, x^2, x^3, x^4]^T](https://s0.wp.com/latex.php?latex=%5Cphi%28x%29+%3D+%5B1%2C+x%2C+x%5E2%2C+x%5E3%2C+x%5E4%5D%5ET&bg=ffffff&fg=444444&s=0&c=20201002)
and the resulting 
would therefore have five columns instead of two (![[1, x]^T](https://s0.wp.com/latex.php?latex=%5B1%2C+x%5D%5ET&bg=ffffff&fg=444444&s=0&c=20201002)
), so now the task is to find ![\boldsymbol{\beta}=[\beta_0, \beta_1, \beta_2, \beta_3, \beta_4]^T](https://s0.wp.com/latex.php?latex=%5Cboldsymbol%7B%5Cbeta%7D%3D%5B%5Cbeta_0%2C+%5Cbeta_1%2C+%5Cbeta_2%2C+%5Cbeta_3%2C+%5Cbeta_4%5D%5ET&bg=ffffff&fg=444444&s=0&c=20201002)
. Following the same logic as before for 
and 
, where prime denotes the new set random variable 
from function 
Water Programming: A Collaborative Research Blog 