Disrupt Spatial Analysis
I have an idea. Hear me out…
If you’re accustomed with this field you probably know it didn’t change that much over the years.
Most geospatial problems involve datasets about geometries, for simplicity sake let’s assume points represented on the surface of the earth, and some metric or query that has to be calculated. Frequently these points will be represented in degrees according to the WGS84 coordinate system (used in GPS). The most common metric to be calculated is Geodesic Distance and the most common query performed is Nearest Neighbor: given two sets of points, A and B, for each point in A get the closest point from B.
For calculating distances there are usually three options:
- Haversine forumla, or great-circle distance, that determines the shortest distance between two points on the surface of a sphere. Even though the earth is not exactly a sphere, this metric allows for a very good approximation that will answer 99% of all use cases. The main advantage is that it is relatively fast to calculate and can be vectorized in NumPy. The main disadvantage is that to get the distance in meters or another unit you will have to multiply the result with the radius of the earth, and frequently different implementations use different radius constants, for example GeoPy uses 6371.009 Km. This means that even though it is consistent it may be not accurate enough for some use cases.
- Vincenty’s formulae, considers the earth to be an ellipsoid, it is extremely accurate, and all it’s parameters are standardized in the WGS84 ellipsoid model. The tradeoff is computing cost, because this method is iterative it is orders of magnitude slower than Haversine. Keep in mind that this method will also results in an approximation, given that the real shape of the earth looks more like a “lumpy potato”, a well placed ellipse still does not match the surface perfectly. For some reason mathematicians have always had a hard time with ellipses, so getting a method that can accurately calculate distances with this model of the surface of the earth was apparently not easy (1975). In Python there’s Geopy or GeographicLib to calculate this distance. GeographicLib in particular offers a 5 nanometer accuracy guarantee.
- Projections, imagine the earth is an orange, you peel the skin of the orange and place it flat on a surface. Now the earth is on a Cartesian plane and distances are euclidean! It is extremely cheap, but comes at a cost of distortion. Projecting a 3d surface in a 2d plane will cause all sorts of problems. There are many projections and some are extremely old, some preserve distances, angles, areas and so on. To address the distortion issue and make projections useful for spatial analysis they are usualy applied to a specific region of earth, providing an accurate representation at the cost of distorting everything that is not on that region. Some examples are ETRS89 for Europe where each country will have it’s own adjusted version or the Universal Transverse Mercator (UTM) that splits the surface of the earth in tiles. In Python we can use PyProj for Coordinate Reference System (CRS) transformations. Libraries like GeoPandas, Shapely and PyGeos require geometries to be in projected coordinate systems in order for most operations to work properly (measurements, predicates, set operations, etc.). Even though it’s possible to calculate distances, angles and areas using the Ellipsoid Model as reference it is faster and easier to do such calculations on a Cartesian plane instead of using non-eucledian geometry.
For Nearest Neighbor queries that’s another story.
- The simplest approach is Brute Force given two sets of points A and B, let’s calculate all possible distances between points and filter the closest pairs. This means there will be len(A) * len(B) total distances to calculate, if the both sets are large the task will be unpractical. The cost of Brute Force approaches is polynomial complexity.
- The solution is using a data structure like a KD-Tree, Ball-Tree or R-Tree to perform the query efficiently. The distance problem still remains, these structures are usually designed to support euclidean or L2 metric. Skleran’s Ball-Tree supports Haversine, however L2 will always be faster.
And then I started thinking of NLP and a embedding model called Skip-gram. In the field of NLP in recent years the most favored representation of a word is an embedding, or a vector of n dimensions, that represents that word in a space. If you calculate the Cosine distance between word vectors it’s possible to find all sorts or relations. So instead of using sparse one-hot-encoding or TF-IDF vectors with high dimensionality you get this beautiful dense and compact vector representation embedded in a space where all words of your vocabulary are magically organized according to their properties and relationship with each other.
Implementing a Skip-gram model in Keras is easy:
- Sample word pairs that are in the same context (positive examples) and word pairs that are not (negative examples);
- Use the embedding layer in Keras so each word can be assigned a vector of random weights;
- Input word pairs to the network, and encode each word form each pair with the same embedding (weight sharing);
- Calculate the cosine distance of both vectors (normalized dot product);
- Optimize embedding weights with Backpropagation and MSE cost;
Why don’t we apply the same principle to the surface of the earth?
Let’s make a projected coordinate system with a Neural Network, here’s the catch let’s use more than 2 dimensions. So we want to project the surface of the earth in a n-dimensional cartesian space, where the l2-norm between two vectors is the Geodesic Distance.
We sacrifice our ability to visualize the points in space, and it requires a slightly more memory, but it should be extremely efficient in terms of computation time, since there are many libraries that support multidimensional Nearest Neighbors with L2 metric like Facebook’s FAISS or SciPy’s KDtree. Unlike traditional map projections we want to make sure that this one works for the whole earth and not only locally and that it is almost as accurate as Vincenty’s method at least from zero meters up to an upper bound.
If this concept is possible the next step would be “hashing” other geometries including Lines and Polygons using this representation has base. Similarly to how sentences are encoded in NLP by merging the embedding vectors of each word. Most operations supported by Shapely and Pygeos would then be implemented and supported.
Imagine having a no hassle, highly accurate and extremely fast and geospatial analysis toolkit mostly based on Neural Networks.
Steps
1 — Sample random pairs of starting points and angles and calculate distances using GeographicLib’s Direct method. It takes as input Latitude and Longitude coordinates of a starting point, a direction in degrees and a distance in meters and outputs the coordinates of a destination point. The main advantage of this approach is that it allows to control the target’s distribution. To speed up the process I used DASK.
2 — Project the points to ECEF (Earth-centered, Earth-fixed coordinate system). Instead of inputing angles into the network we are going to input Cartesian coordinates of the surface of the earth, in a plane in which the origin is the center of the earth. This step is extremely cheap.
Another advantage of ECEF is that it already does some of the work for us, since the surface has low curvature for points closer together this means we can calculate accurate distances up to an upper bound (e.g. ≤10-100 KM).
3—Create the Neural Network with Keras.
The basic logic here is similar to an embedding model. The input to the network is pairs of points (n, 2, 3) and the output is the distance in meters. We have a “base model” that can have whatever layers we want, we apply this base model to the point on the left and to the point on the right (weight sharing), then the output representation a vector of 3 dimensions is concatenated with the original ECEF coordinates. This means each projected coordinate will have 6 dimensions. After that we simply calculate the L2 Norm of both vectors, implemented using a custom Keras Lambda layer, getting the predicted distance. We then optimize the network with MSE trying to reduce error between predicted and true distance.
Because of numerical instability issues the L2 Norm in Tensorflow has to be written with a custom gradient.
It is very important to set Keras precision to 64 bit floats, otherwise there won’t be enough precision to represent ECEF coordinates, consequently small distances will be misrepresented.
Keep in mind that the scale of the error has to be proportional with distance, this means we want to achieve a model that has very low error on low distances and also low error on long distances.
The code for the network is shared below. If someone wants the sampling notebook please request in the comments and I will publish it. My sample was approximately 60 million random points on the surface of the earth, distances range from 0 meters to 1 million meters.
Is this working?
No, I was able to achieve some promising results but there’s always some kind of drawback. I’m sharing this hoping someone can prove or disprove the idea.
Proof of Concept: https://carlostgameiro.medium.com/projections-using-neural-networks-a812e98777d
