Introducing kNN Join for Wherobots and Apache Sedona Posted on November 5, 2024October 4, 2026 by Ben Pruden Editor’s note: The Wherobots Spatial Data Catalog is now the Havasu Catalog. TL;DR: WherobotsDB and Apache Sedona 1.7.0 introduce two types of kNN Join for large-scale geospatial workloads. The Exact kNN Join (ST_KNN) returns precise nearest neighbors and is best when accuracy is critical. The Approximate kNN Join (ST_AKNN) trades roughly 3-7% accuracy for significantly faster query times and lower cost. Both outperform PostGIS at scale, with PostGIS unable to complete joins at the 10M x 1B dataset size. We are excited to introduce the k-Nearest Neighbors Join (kNN Join) for WherobotsDB and Apache Sedona 1.7.0. With the kNN Join, you can efficiently find the closest entities to your points of interest from datasets at scale. We’ve released two types of kNN Joins: the Exact kNN Join and the Approximate kNN Join. With the Exact kNN Join you get precise results, but the additional precision adds cost and time to the query. The Exact kNN Join is useful when precise object proximity is critical. It is now available in WherobotsDB and will be available in Apache Sedona 1.7.0. The Approximate kNN Join delivers results faster and at a lower cost than the Exact kNN Join, and is tested to be ~3-7% less accurate than the Exact kNN Join for varying data sizes (more info below). The Approximate kNN Join, available now in WherobotsDB, is useful when performance matters more than accuracy, such as testing workloads or when query speed matters more than precision. Click here to launch this interactive notebook Launch Notebook When is the kNN Join useful? The kNN Join allows you to quickly find a chosen set (k) of nearest entities (objects) to your points of interest (queries). When Should You Use the Exact kNN Join The Exact kNN Join will efficiently identify the exact k closest objects to queries of interest. This makes it a highly useful algorithm when precise object proximity is crucial. Common applications of the Exact kNN Join at scale include: Reverse geocoding: Ride share services can use the Exact kNN Join to map high quantities of user tracking coordinates to street address for smooth rider pickups. Identifying nearest transit options: Map services can use the Exact kNN Join to accurately locate (or simulate) the closest public transit options for many people simultaneously. Analyzing the distribution of public services: City planners can use the Exact kNN Join to identify geographic areas without public services nearby to inform future development plans. Efficient freight routing: Shipping and cargo trucks can use the Exact kNN Join to locate exact locations of nearest refueling stations, ramps, marine life, etc. to create efficient routes and reduce time and cost. When Should You Use the Approximate kNN Join The Approximate kNN Join trades off accuracy for higher query speed. It will find a set of k close neighbors (but not necessarily the closest neighbors) to your query set. Some common applications of the Approximate kNN Join include: Local search: Map services use Approximate kNN Join to help users find places of interest within a certain distance of their geolocation. Image search: Geospatial modeling teams leverage Approximate kNN to find images for geographic points of interest to augment available data for model training. More importantly, the Approximate kNN Join unlocks kNN to be used for any compute or time constrained workloads, including: Testing: Teams developing geospatial workloads leverage the speed of the Approximate kNN Join to quickly test and iterate their pipelines. Meeting latency requirements: Teams use Approximate kNN Join to meet latency requirements and serve identified nearest neighbors at scale in production. Why Self-Hosting kNN Join Algorithms Is Harder Than It Looks Before, teams leveraging the kNN algorithm for geospatial workloads would have to: Determine how to host and scale general-purpose open source kNN solutions, including managing the curse of dimensionality, manage computational complexity that increases with data size, and preprocessing outlier data. Limit the frequency of executing kNN join workloads due to high time and compute costs. What You Gain by Using kNN Join in Wherobots or Apache Sedona Now with the kNN Join in Wherobots or Apache Sedona, you can: Leverage geospatially optimized kNN joins that scale efficiently with increasing data sizes without the overhead of building and maintaining your own solution. Iterate at a higher pace, by leveraging the Approximate kNN Join to trade-off performance versus accuracy. And, as your workload requires, run the Exact kNN Join to guarantee accurate neighbors for points of interest. See how ST_KNN and ST_AKNN perform at 44 million geometries across five US states. We’ll walk through a brief overview of the kNN Join, share how each perform at various scale, and how to use it. How kNN Join Works in WherobotsDB and Apache Sedona When used with geospatial data, the kNN Join identifies the k-nearest neighbors for a given spatial point or region based on geographic proximity. It involves two datasets: a query dataset containing points of interest, and an objects dataset, containing the entities you’re trying to locate near those points of interest. For each record in the queries dataset, the kNN Join finds the resulting set (k) of nearest neighbors from the objects dataset based on a user-defined distance metric and k value. For example, let’s say you’re interested in locating the three closest coffee shops to specific public transit stops in your neighborhood. The query dataset would be gps coordinates of the public transit stops. The objects dataset would be gps coordinates of all coffee shops in your neighborhood. As we’re looking for the three closest coffee shops, the k value is 3. We’ll say that our distance metric is three blocks. Using this as input, the kNN Join will conduct a search of the objects and return three coffee shops that are within three blocks of each of the public transit stops in your neighborhood.An illustration of kNN join. Using the example above, the red dots are the query dataset containing gps coordinates of the public transit stops. The green dots are the objects dataset containing gps coordinates of all coffee shops in your area. We want the 3-nearest coffee shops (k=3) to these public transit stops. Using kNN join, we find the closest coffee shops to the public transit stops. The query dataset and only their closest found neighbors are shown in the df_joined. kNN Join SQL API Overview We’ve made it easy to integrate the kNN Join into your workflows. Here is a quick overview of the kNN Join Spatial SQL APIs. For an in-depth look, see our docs here. Exact kNN Join Supported GeometriesHyperparametersOutputSQL APIExact kNN Joinpoints, geometriesnumber of nearest neighbors to identify (k)dataframe with k nearest neighbors per query ranked by distanceST_KNN(...) Exact_kNN_Join_Results = ST_KNN(WEATHER.GEOMETRY, FLIGHTS.GEOMETRY, 4, FALSE) Approximate kNN Join Supported GeometriesHyperparametersOutputSQL APIApproximate kNN Joinpoints, geometriesnumber of nearest neighbors to identify (k)dataframe with k nearest neighbors per query ranked by distanceST_AKNN(...) Approximate_kNN_Join_Results = ST_AKNN(WEATHER.GEOMETRY, FLIGHTS.GEOMETRY, 4, FALSE) kNN Join Performance Benchmarks: Wherobots vs PostGIS We’ll compare how the Exact and Approximate kNN Joins perform vs PostGIS, and how the Approximate kNN Join performs against the Exact kNN Join. kNN Join vs PostGIS Let’s look at some performance benchmarks for the Exact and Approximate kNN Joins compared to PostGIS’s exact kNN join algorithm. The following benchmarks are based on measured wall clock time for each join on varying data sizes. Each data size is defined by the query data size and the object data size, so a query set of 1M records and a object set of 12M records is described as 1M x 12M in the graph below. To compare scale and performance, we avoid compute comparisons across options by normalizing on cost per hour for each join. As seen below, the Approximate and Exact kNN Joins scale better than PostGIS with increasing dataset sizes. A key limitation of PostGIS’s exact kNN join is its inability to scale horizontally across machines, which is why performance with large scale exact kNN joins is better with Wherobots. Meanwhile, the Exact and Approximate kNN Joins distribute the join horizontally across workers and process larger datasets more efficiently. Specifically, the Approximate kNN Join is more performant and efficient at cost per workload than PostGIS and Exact kNN Join. Meanwhile, the Exact kNN Join is more effective for larger workloads than PostGIS, with PostGIS unable to complete the join for the 10M x 1B data size. Approximate vs Exact KNN Join Benchmarks Let’s look at some performance benchmarks to help you choose which kNN Join is right for a given workload. We’ve chosen a variety of metrics to evaluate the accuracy of the Approximate kNN Join. Accuracy: Measures the percentage of the Approximate kNN Join results that are the closest objects to queries. Jaccard Index: Measures the overlap between the Approximate kNN Join results and the Exact kNN Join results. Let’s take a look at wall clock time for the Approximate kNN Join versus Exact kNN Join for increasing data sizes. Let’s take a look at accuracy and Jaccard scores for the Approximate kNN Join for increasing data sizes. As seen above, the Approximate kNN Join will run faster than the Exact kNN Join for datasets of varying size, is highly accurate, and identifies highly similar sets of nearest neighbors as the Exact KNN Join, even for increasing dataset sizes. We recommend leveraging the Approximate kNN Join whenever you’re workloads allow for less precision, are cost constrained, or require high performance. Tutorial: How to Run kNN Join on Weather and Flight Data in WherobotsDB The following example details how to apply the Exact and Approximate kNN Joins on weather events data to find which flights are closest to each weather event, which can be crucial for making real-time decisions in air traffic management and ensuring flight safety. We will load this weather data from the Wherobots Spatial Catalog. The Spatial Catalog includes a collection of Wherobots maintained open datasets from various data sources like Overture Maps, LandSAT, Wild Fires, New York Taxi Data, and more. These datasets are optimized for fast and efficient analytics with WherobotsDB. The queries table contains the locations of weather events, such as storms or turbulence, while the objects table contains flight locations. Visit the user documentation for a full walk through of both kNN joins. 1. Preparing the Weather Events Data as Queries DataFrame First, we load the weather events data from the Spatial Data Catalog. Then, we assign a monotonically increasing ID to each row in the weather events, and load it to a temporary view so we can use in the SQL. df_queries = sedona.table("wherobots_pro_data.weather.weather_events") df_queries = df_queries.withColumn("id", monotonically_increasing_id()) df_queries.createOrReplaceTempView("weather") df_queries.select("id", "geometry", "type").show(20, False) The weather queries events DataFrame looks like this: +-------+-------------------------+----+ |id |geometry |type| +-------+-------------------------+----+ |2716677|POINT (-110.974 41.8191) |Snow| |2656934|POINT (-109.4604 41.5323)|Snow| |2653971|POINT (-109.0528 41.5945)|Rain| |2519252|POINT (-105.0333 42.45) |Rain| |2480600|POINT (-105.5419 44.3394)|Snow| +-------+-------------------------+----+ only showing top 5 rows 2. Preparing the Flights Data as Objects DataFrame We are using the flights tracking data from ABS-B flight observation data. This data originated from ADSB.lol, which provides publicly available daily updates of crowdsourced flight tracking data. df_objects = sedona.read.format("geoparquet").load("s3a://wherobots-examples/data/examples/flights/2024_s2.parquet") df_objects.createOrReplaceTempView("flights") The flights objects DataFrame looks like this: +---------------------------+-----------------------------+-------------------+ |desc |geometry |timestamp | +---------------------------+-----------------------------+-------------------+ |BOEING 737-800 |POINT (-110.232108 43.929066)|2024-01-11 00:00:00| |CESSNA 240 CORVALIS TTX |POINT (-110.852247 43.765574)|2024-01-01 00:00:00| |CESSNA 172 SKYHAWK |POINT (-108.790218 44.801926)|2024-01-23 00:00:00| |CESSNA 510 CITATION MUSTANG|POINT (-110.795842 43.518448)|2024-01-26 00:00:00| |BOEING 737-800 |POINT (-110.873321 43.079086)|2024-01-18 00:00:00| +---------------------------+-----------------------------+-------------------+ only showing top 5 rows 3. Run the kNN Join Once the two datasets are loaded, we can use the following SQL to run the kNN Joins. The query below runs the Exact kNN Join with the number of neighbors set to 4. df_knn_join = sedona.sql(""" SELECT WEATHER.GEOMETRY AS WEATHER_GEOM, FLIGHTS.GEOMETRY AS FLIGHTS_GEOM FROM WEATHER JOIN FLIGHTS ON ST_KNN(WEATHER.GEOMETRY, FLIGHTS.GEOMETRY, 4, FALSE) """) You can also run the Approximate kNN Join with number of neighbors set to 4, and add UDFs to the result DataFrame as following. df_knn_join = sedona.sql(""" SELECT WEATHER.GEOMETRY AS WEATHER_GEOM, [WEATHER.ID](http://weather.id/) AS QID, FLIGHTS.GEOMETRY AS FLIGHTS_GEOM, ST_DISTANCESPHERE(WEATHER.GEOMETRY, FLIGHTS.GEOMETRY) AS DISTANCE, ST_MAKELINE(WEATHER.GEOMETRY, FLIGHTS.GEOMETRY) AS LINE FROM WEATHER JOIN FLIGHTS ON ST_AKNN(WEATHER.GEOMETRY, FLIGHTS.GEOMETRY, 4, FALSE) """) The joined DataFrame looks like the following for a single query point. +-------+-------------------------+-----------------------------+-----------------+ |QID |QUERIES_GEOM |OBJECTS_GEOM |DISTANCE | +-------+-------------------------+-----------------------------+-----------------+ |2804966|POINT (-110.4211 44.5444)|POINT (-110.425023 44.544205)|311.6515627672718| |2804966|POINT (-110.4211 44.5444)|POINT (-110.416483 44.546535)|436.1578754802839| |2804966|POINT (-110.4211 44.5444)|POINT (-110.428809 44.546302)|646.4967456258582| |2804966|POINT (-110.4211 44.5444)|POINT (-110.415679 44.550445)|797.7249421529948| +-------+-------------------------+-----------------------------+-----------------+ only showing top 4 rows+-------+-------------------------+-----------------------------+-----------------+ |QID |QUERIES_GEOM |OBJECTS_GEOM |DISTANCE | +-------+-------------------------+-----------------------------+-----------------+ |2804966|POINT (-110.4211 44.5444)|POINT (-110.425023 44.544205)|311.6515627672718| |2804966|POINT (-110.4211 44.5444)|POINT (-110.416483 44.546535)|436.1578754802839| |2804966|POINT (-110.4211 44.5444)|POINT (-110.428809 44.546302)|646.4967456258582| |2804966|POINT (-110.4211 44.5444)|POINT (-110.415679 44.550445)|797.7249421529948| +-------+-------------------------+-----------------------------+-----------------+ only showing top 4 rows 4. Visualize the kNN Join Results We can also visualize the Exact kNN Join results using SedonaKepler by adding the joined DataFrame to a map. # create map for the results map_view = SedonaKepler.create_map(df_unique_qid.select('QUERIES_GEOM'), name="WEATHER EVENTS") SedonaKepler.add_df(map_view, df=df_related_rows.select('OBJECTS_GEOM', 'DISTANCE').withColumnRenamed("OBJECTS_GEOM", "geometry"), name="FLIGHTS") SedonaKepler.add_df(map_view, df=df_related_rows.select('LINE', 'DISTANCE').withColumnRenamed("LINE", "geometry"), name="KNN LINES") # show the map map_view Below we see a portion of the map we generated above. In this portion of the map, we see the results of the Exact kNN Join. Given our distance requirement between a weather event and flight (the lines), for each of the weather events from our query database (blue dots), we see the 4 closest flights to each weather event (yellow dots) . Get started with the Wherobots kNN Join Wherobots Users All Wherobots users have access to the Exact and Approximate kNN Joins. If you haven’t already, create a free Wherobots organization subscribed to the Community Edition of Wherobots. Start a Wherobots Notebook In the Notebook environment, explore the python/wherobots-db/wherobots-db-knn-joins.ipynb example that you can use to get started. Need additional help? Check out our user documentation, and send us a note if needed at support@wherobots.com. What’s next We’re excited to hear what algorithms you’d like us to support. We can’t wait for your feedback and to see what you’ll create! Start building with Wherobots Get Started Key takeawaysWherobotsDB and Apache Sedona 1.7.0 add two k-nearest-neighbor joins: Exact kNN (ST_KNN) for precise neighbors, and Approximate kNN (ST_AKNN) when speed and cost matter more.ST_AKNN is tested at about 3–7% less accurate than ST_KNN across varying data sizes, with Jaccard overlap also reported against the exact results.Both scale better than PostGIS kNN because they distribute horizontally. PostGIS could not complete the join at the 10 million × 1 billion dataset size in the benchmark.The post points to a ST_KNN / ST_AKNN example on 44 million geometries across five US states, and a tutorial joining weather events to ADS-B flights with k=4.Exact and Approximate kNN Joins are available to all Wherobots users, including Community Edition, via the python/wherobots-db/wherobots-db-knn-joins.ipynb notebook.
Introducing GeoStats for WherobotsAI and Apache Sedona Posted on October 8, 2024October 3, 2026 by Ben Pruden Introducing GeoStats for WherobotsAI and Apache Sedona We are excited to introduce GeoStats, a machine learning (ML) and statistical toolbox for WherobotsAI and Apache Sedona users. With GeoStats, you can easily identify critical patterns in geospatial datasets such as hotspots and anomalies, and quickly get critical insights from large scale data. While these algorithms are supported in other packages, we’ve optimized each algorithm to be highly performant for small to planetary scale geospatial workloads. That means, you can get results from these algorithms significantly faster, at a lower cost, and do it all more productively, through a unified development experience purpose-built for geospatial data science and ETL. The Wherobots toolbox supports DBSCAN, Local Outlier Factor (LOF), and Getis-Ord (Gi*) algorithms. Apache Sedona users can utilize DBSCAN starting with Apache Sedona version 1.7.0, and like all other features of Apache Sedona, its fully compatible with Wherobots. Use Cases for GeoStats DBSCAN DBSCAN is the most popular algorithm we see in geospatial use cases. It identifies clusters, areas of your data that are closely packed together, and outliers, areas of your data that are set apart. Typical use cases for DBSCAN are found in: Retail: Decision makers use DBSCAN with location data to understand areas of high and low pedestrian activity to decide where to setup retailing establishments. City Planning: City planners use DBSCAN with GPS data to optimize transit support by identifying high usage routes, areas in need of additional transit options, or areas that have too much support. Air Traffic Control: Traffic controllers use DBSCAN to identify areas with increasing weather activity to optimize flight routing. Risk computation: Insurers and others can use DBSCAN to make policy decisions and calculate risk where risk is correlated to the proximity of two or more features of interest. Local Outlier Factor (LOF) LOF is an anomaly detection algorithm that identifies outliers present in a dataset. Typical use cases for LOF include: Data analysis and cleansing: Data teams can use LOF to identify and remove anomalies within a dataset, like removing erroneous GPS data points from a trace dataset Getis-Ord (Gi*) Getis-Ord is also a popular algorithm for identifying local hot and cold spots. Typical use cases for Gi* include: Public health: Officials can use disease data with Gi* to identify areas of abnormal disease outbreak Telecommunications: Network administrators can use Gi* to identify areas of high demand and optimize network deployment Insurance: Insurers can identify areas prone to specific claims to better manage risk Click here to launch this interactive notebook Launch Notebook Traditional challenges with using these algorithms on geospatial data Before GeoStats, teams leveraging any of the algorithms in the toolbox in a data analysis or ML pipeline would: Struggle to get performance or scale from the underlying solutions that also don’t perform well when joining geospatial data. Determine how to host and scale open source versions of popular ML and statistical algorithms, like PostGIS or scikit-learn DBSCAN, PySal Gi*, or scikit-learn LOF, to work for geospatial data types and geospatial data formats. Replicate this overhead each time they want to deploy a new algorithm for geospatial data. Benefits of WherobotsAI GeoStats With GeoStats in WherobotsAI, you can now: Easily run native algorithms on a cloud-based engine, optimized for producing spatial data products and insights at scale. Use these algorithms without the operational overhead associated with setup and maintenance. Leverage optimized, hosted algorithms within a single platform to easily experiment and get critical insights faster. We’ll walk through a brief overview of each algorithm, how to use them, and show how they perform at various scales. Diving Deeper into the GeoStats Toolbox DBSCAN Overview DBSCAN is a density-based clustering algorithm. Given a set of points in some space, it groups points with many nearby neighbors and marks as outlier points that lie alone in low-density regions. How to use DBSCAN in Wherobots The following examples assume you have already setup an organization and have an active runtime and notebook, with a dataframe of interest to run the algorithms on. WherobotsAI GeoStats DBSCAN Python API OverviewFor a full walk through see the Python API reference: dbscan(...). Supported Geometries: points, linestrings, polygons Hyperparameters: max distance to neighbors (epsilon), min neighbors (min_points) Output: dataframe with cluster id DBSCAN Walk Through Choose your dataset and create a Sedona DataFrame. dataset=sedona.createDataFrame(X).select(ST_MakePoint("_1", "_2").alias("geometry")) Choose values for your hyperparameters, max distance to neighbors (epsilon) and minimum neighbors (min_points). These values will determine how DBSCAN identifies clusters. epsilon=0.3 min_points=10 Run DBSCAN on your DataFrame with your chosen hyperparameter values. clusters_df = dbscan(df, epsilon=0.3, min_points=10, include_outliers=True) Analyze the results. For each datapoint, DBSCAN returns the cluster it’s associated with or if it’s an outlier. +--------------------+------+-------+ | geometry|isCore|cluster| +--------------------+------+-------+ |POINT (1.22185277...| false| 1| |POINT (0.77885034...| false| 1| |POINT (-2.2744742...| false| 2| +--------------------+------+-------+ only showing top 3 rows There’s a complete example of how to use DBSCAN in the Wherobots user documentation. DBSCAN Performance Overview To show DBSCAN performance in Wherobots, we created a European sample of the Overture buildings dataset, and ran DBSCAN to identify clusters of buildings near each other, starting from the geographic center of Europe and worked outwards. For each subsampled dataset, we run DBSCAN with an epsilon of 0.005 degrees (i.e. ~30 feet) and min_points value of 4 on a Large runtime in Wherobots Cloud. As seen below, DBSCAN effectively processes an increasing number of records, with 100M records taking 1.6 hrs to process. Local Outlier Factor (LOF) LOF is an anomaly detection algorithms that identifies outliers present in a dataset. It does this by measuring how close a given data point is to a set of k-nearest neighbors (with k being a user chosen hyperparameter) in comparison to how close its nearest neighbors are to their nearest neighbors. LOF provides a score that represents the degree to which a record is an inlier or outlier. How to use LOF in Wherobots For the full example, please see this docs page. WherobotsAI GeoStats LOF Python API OverviewFor a full walk through see the Python API reference: local_outlier_factor(...). Supported Geometries: points, linestrings, polygons Hyperparameters: number of nearest neighbors to use Output: score representing degree of inlier or outlier LOF Walk Through Choose your dataset and create a Sedona DataFrame. df = sedona.createDataFrame(X).select(ST_MakePoint(f.col("_1"), f.col("_2")).alias("geometry")) Choose your k value for how many nearest neighbors you want to use to measure density near a given datapoint. k=20 Run LOF on your DataFrame with your chosen k value. outliers_df = local_outlier_factor(df, k=20) Analyze your results. LOF returns a score for each datapoint representing the degree of inlier or outlier. +--------------------+------------------+ | geometry| lof| +--------------------+------------------+ |POINT (-1.9169927...|0.9991534865548664| |POINT (-1.7562422...|1.1370318880088373| |POINT (-2.0107478...|1.1533763384772193| +--------------------+------------------+ only showing top 3 rows There’s a complete example of how to use LOF in the Wherobots user documentation. LOF Performance Overview We followed the same procedure with DBSCAN but ran LOF to identify clusters of buildings near each other. With each set of buildings we ran LOF with a k=20 on a large Wherobots Cloud runtime. As seen below, GeoStats LOF scales effectively with growing data size with 100M records taking 10 mins to process. Getis-Ord (Gi*) Overview Getis-Ord is an algorithm for identifying statistically significant local hot and cold spots. How to use GeoStats Gi* WherobotsAI GeoStats Gi* Python API OverviewFor the full example, please see this docs g_local(...). Supported Geometries: points, linestrings, polygons Hyperparameters: star, neighbor weighting Output: Set of statistics that indicate the degree of local hot or cold spot for a given record Gi* Walk Through Choose your dataset and create a Sedona Dataframe. places_df = ( sedona.table("wherobots_open_data.overture_2024_07_22.places_place") .select(f.col("geometry"), f.col("categories")) .withColumn("h3Cell", ST_H3CellIDs(f.col("geometry"), h3_zoom_level, False)[0]) ) Choose how you’d like to weight datapoints (ex: datapoints in a specific geographic area need to be weighted higher or any datapoint close to a given datapoint need to be weighted higher) and star (boolean to indicate if a record is a neighbor of itself). star = True neighbor_search_radius_degrees = 1.0 variable_column = "myNumericColumnName" weighted_dataframe = add_binary_distance_band_column( df, neighbor_search_radius_degrees, include_self=star ) Run Gi* on your DataFrame with your chosen hyperparameters. gi_df = g_local( weighted_dataframe, variable_column, star=star ) Analyze your results. For each datapoint, Gi* returns a set of statistics that indicate the degree of local hot or cold spot. +----------+-------------------+--------------------+--------------------+------------------+--------------------+ |num_places| G| EG| VG| Z| P| +----------+-------------------+--------------------+--------------------+------------------+--------------------+ | 871| 0.1397485091609774|0.013219284603421462|5.542296862370928E-5|16.995969941572465| 0.0| | 908|0.16097739240211956|0.013219284603421462|5.542296862370928E-5|19.847528249317246| 0.0| | 218|0.11812096144582315|0.013219284603421462|5.542296862370928E-5|14.090861243071908| 0.0| +----------+-------------------+--------------------+--------------------+------------------+--------------------+ only showing top 3 rows There’s a complete example of how to use Gi* in the Wherobots user documentation. Getis-Ord Performance Overview To showcase how Gi performs in Wherobots, again we used the same example as DBSCAN, but ran Gi on the area of the buildings. With each set of buildings we ran Gi* with a binary neighbor weight and a neighborhood radius of .007 degrees (~0.5 miles) on a Large runtime in Wherobots Cloud. As seen below, the algorithm scales mostly linearly with the number of records, with 100M records taking 1.6 hours to process. Get started with WherobotsAI GeoStats The way we implemented these algorithms for large scale geospatial workloads, will help you make sense of your geospatial data faster. You can get started for free today. If you haven’t already, create a free Wherobots Organization subscribed to the Community Edition of Wherobots. Start a Wherobots Notebook In the Notebook environment, explore the notebook_example/python/wherobots-ai/ folder for examples that you can use to get started. Need additional help? Check out our user documentation, and send us a note if needed at support@wherobots.com. Apache Sedona Users Apache Sedona users will have access to GeoStats DBSCAN with the Apache Sedona 1.7.0 release. Subscribe to the Spatial Intelligence Newsletter and join the Sedona community to get notified of the release and get started! What’s next We’re excited to hear what ML and statistical algorithms you’d like us to support. We can’t wait for your feedback and to see what you’ll create! Start building with Wherobots Access Here Key takeawaysGeoStats is a machine-learning and statistics toolbox for WherobotsAI and Apache Sedona, with DBSCAN clustering, Local Outlier Factor (LOF) anomaly detection, and Getis-Ord Gi* hot/cold spots.Apache Sedona users get DBSCAN starting in Sedona 1.7.0; the Wherobots toolbox includes all three algorithms and stays Sedona-compatible.On a Large Wherobots runtime, DBSCAN on a growing European Overture buildings sample (epsilon 0.005 degrees, about 30 feet, min_points 4) processed 100 million records in 1.6 hours.LOF on the same buildings progression with k=20 processed 100 million records in 10 minutes on a Large runtime.Gi* with binary neighbor weights and a 0.007-degree (about 0.5 mile) radius processed 100 million records in 1.6 hours on a Large runtime, scaling mostly linearly.
Announcing SAML Single Sign-On (SSO) Support Posted on September 25, 2024October 3, 2026 by Pranav Toggi We’re excited to announce that Wherobots Cloud now supports SAML SSO for Professional Edition customers. This enhancement underscores our commitment to providing robust security features tailored to the needs of businesses that handle sensitive data. What is SAML Single Sign-On Many companies already use SAML SSO. It is the mechanism by which you log in to third party tools with your company’s centralized login (Google, Okta, etc.) SAML (Security Assertion Markup Language) are essential components of modern security protocols. SAML is an open standard for exchanging authentication and authorization data between parties and SSO (Single Sign-On) allows users to log in once and gain access to multiple systems without re-entering credentials. How to enable SAML SSO Setting up SAML SSO is straightforward. Any Professional Edition Wherobots administrator can enable this feature by following these steps: Verify Your Domain:Go to your organization’s settings, copy the provided TXT DNS record, and add it to your domain provider’s DNS settings. Once done, click “Verify” in Wherobots. Configure Your Identity Provider (IdP):Log in to your IdP (e.g., Google Workspace, Okta, OneLogin, Azure AD) and configure it using the SAML details provided in the Wherobots settings. Enter IdP Details into Wherobots:Input the Identity Provider Issuer URL, SSO URL, and certificate details from your IdP into the SAML section in Wherobots. Enable SAML SSO:Make sure there’s an admin user with an email from the verified domain. Then, toggle the “Enable” switch in Wherobots to activate SSO. Test the Integration:Log in using your verified domain email to ensure everything is set up correctly. For more detailed step-by-step instructions, check out our comprehensive getting started guide. Why Your Organization Should Use SAML Single Sign-On Benefits to Users For users, SAML SSO simplifies the login process. With fewer passwords to remember and less frequent logins, users save time and reduce frustration. This streamlined access means employees can focus more on their tasks and less on managing multiple credentials. Benefits to Organizations Organizations benefit from SAML SSO by enhancing security and reducing administrative overhead. Centralized authentication means fewer password-related support tickets and a lower risk of phishing attacks, as users are less likely to fall for credential-stealing schemes. Moreover, it ensures compliance with security policies and regulatory requirements by enforcing strong, consistent access controls. It also allows organizations to grant access to certain third party services only for certain users/groups of users. Location data is particularly sensitive, as it can reveal personal habits, routines, and even confidential business operations. For example, in healthcare, location data of patients visiting specialized clinics could inadvertently disclose medical information. By implementing SAML SSO, organizations can better control access, ensuring that only authorized personnel can view and interact with this information within the Wherobots platform. Elevate Your Security with SAML SSO Take this opportunity to simplify authentication, reduce security risks, and improve productivity. If you are not already a Professional Edition customer, upgrade to gain access to SAML SSO. By upgrading, you’ll not only bolster your security measures and simplify logins, but you’ll gain access to more powerful workloads, service principals, and more. Contact us to schedule a demo or get started! Want to keep up with the latest developer news from the Wherobots and Apache Sedona community? Sign up for the Spatial Intelligence Newsletter: Key takeawaysWherobots Cloud now supports SAML SSO for Professional Edition customers, aimed at organizations that already centralize login through an identity provider.SAML exchanges authentication and authorization data between parties; SSO lets users log in once and reach multiple systems without re-entering credentials.Any Professional Edition administrator can enable it: verify the company domain with a TXT DNS record, configure the IdP, paste issuer URL, SSO URL, and certificate into Wherobots, ensure an admin email on the verified domain, then toggle Enable.Documented identity providers include Google Workspace, Okta, OneLogin, and Azure AD.The post argues SSO cuts password-reset tickets and phishing risk and helps control who can see sensitive location data—for example healthcare clinic visits—inside Wherobots. Upgrading to Professional also unlocks more powerful workloads and service principals.
🌶 Comparing taco chains :: a consumer retail cannibalization study with isochrones Posted on June 11, 2024October 4, 2026 by Ben Pruden Authors Ilya Marchenko Using Wherobots for a Retail Cannibalization Study Comparing Two Leading Taco Chains In this post, we explore how to implement a workflow from the commercial real estate (CRE) space using Wherobots Cloud. This workflow is commonly known as a cannibalization study, and we will be using WherobotsDB, POI data from OvertureMaps, the open source Valhalla API, and visualization capabilities offered by SedonaKepler. What is a retail cannibalization study? In CRE (consumer real estate), stakeholders are often interested in questions like “If we build a new fast food restaurant here, how will its performance be affected by other similar fast food locations that already exist nearby?”. The idea of the new fast food restaurant “eating into” the sales of other fast food restaurants that already exist nearby is what is known as ‘cannibalization’. The main objective of studying this phenomenon is to determine the extent to which a new store might divert sales from existing stores owned by the same company or brand and evaluate the overall impact on the company’s market share and profitability in the area. Cannibalization Study in Wherobots For this case study, we will look at two taco chains which are located primarily in Texas: Torchy’s Tacos and Velvet Taco. In general, information about the performance of individual locations and customer demographics are often proprietary information. We can, however, still learn a great deal about the potential for cannibalization both between these two chains as competitors, and between individual locations of each chain. We also know, based on our own experience, these chains compete with each other. Which taco shop to go to when we are visiting Texas is always a spicy debate. We begin by importing modules that will be useful to us as we go on. import geopandas as gpd import pandas as pd import requests from sedona.spark import * from pyspark.sql.functions import explode, array from pyspark.sql import functions as F Next, we can initiate a Sedona context. config = SedonaContext.builder().getOrCreate() sedona = SedonaContext.create(config) Identifying Points of Interest Now, we need to retrieve the locations of Torchy’s Tacos and Velvet Taco locations. In general, one can do this via a variety of both free and paid means. We will look at a simple, free approach that is made possible by the integration of Overture Maps data into the Wherobots environment: sedona.table("wherobots_open_data.overture.places_place"). \\ createOrReplaceTempView("places") We create a view of the Overture Maps places database, which contains information on points of interest (POI’s) worldwide. Now, we can select the POI’s which are relevant to this exercise: stores = sedona.sql(""" SELECT id, names.common[0].value as name, ST_X(geometry) as long, ST_Y(geometry) as lat, geometry, CASE WHEN names.common[0].value LIKE "%Torchy's Tacos%" THEN "Torchy's Tacos" ELSE 'Velvet Taco' END AS chain FROM places WHERE addresses[0].region = 'TX' AND (names.common[0].value LIKE "%Torchy's Tacos%" OR names.common[0].value LIKE '%Velvet Taco%') """) Calling stores.show() gives us a look at the spark DataFrame we created: +--------------------+--------------+-----------+----------+ | id| name| long| lat| +--------------------+--------------+-----------+----------+ |tmp_8104A79216254...|Torchy's Tacos| -98.59689| 29.60891| |tmp_D17CA8BD72325...|Torchy's Tacos| -97.74175| 30.29368| |tmp_F497329382C10...| Velvet Taco| -95.48866| 30.18314| |tmp_9B40A1BF3237E...|Torchy's Tacos| -96.805853| 32.909982| |tmp_38210E5EC047B...|Torchy's Tacos| -96.68755| 33.10118| |tmp_DF0C5DF6CA549...|Torchy's Tacos| -97.75159| 30.24542| |tmp_BE38CAC8D46CF...|Torchy's Tacos| -97.80877| 30.52676| |tmp_44390C4117BEA...|Torchy's Tacos| -97.82594| 30.4547| |tmp_8032605AA5BDC...| Velvet Taco| -96.469695| 32.898634| |tmp_0A2AA67757F42...|Torchy's Tacos| -96.44858| 32.90856| |tmp_643821EB9C104...|Torchy's Tacos| -97.11933| 32.94021| |tmp_0042962D27E06...| Velvet Taco|-95.3905374|29.7444214| |tmp_8D0E2246C3F36...|Torchy's Tacos| -97.15952| 33.22987| |tmp_CB939610BC175...|Torchy's Tacos| -95.62067| 29.60098| |tmp_54C9A79320840...|Torchy's Tacos| -97.75604| 30.37091| |tmp_96D7B4FBCB327...|Torchy's Tacos| -98.49816| 29.60937| |tmp_1BB732F35314D...| Velvet Taco| -95.41044| 29.804| |tmp_55787B14975DD...| Velvet Taco|-96.7173913|32.9758554| |tmp_7DC02C9CC1FAA...|Torchy's Tacos| -95.29544| 32.30361| |tmp_1987B31B9E24D...| Velvet Taco| -95.41006| 29.770256| +--------------------+--------------+-----------+----------+ only showing top 20 rows We’ve retrieved the latitude and longitude of our locations, as well as the name of the chain each location belongs to. We used the CASE WHEN statement in our query in order to simplify the location names. This way, we can easily select all the stores from the Torchy’s Tacos chain, for example, and not have to worry about individual locations being called things like “Torchy’s Tacos – Rice Village” or “Velvet Taco Midtown”, etc. We can also visualize these locations using SedonaKepler. First, we can create the map using the following snippet: location_map = SedonaKepler.create_map(stores, "Locations", config = location_map_cfg) Then, we can display the results by simply calling location_map in the notebook. For convenience, we included the location_map_cfg Python dict in our notebook, which stores the settings necessary for the map to be created with the locations color-coded by chain. If we wish to make modifications to the map and save the new configuration for later use, we can do so by calling location_map.config and saving the result either as a cell in our notebook or in a separate location_map_cfg.py file. Generating Isochrones Now, for each of these locations, we can generate a polygon known as an isochrone or drivetime. These polygons will represent the areas that are within a certain time’s drive from the given location. We will generate these drivetimes using the Valhalla isochrone api: def get_isochrone(lat, lng, costing, time_steps, name, location_id): url = "<https://valhalla1.openstreetmap.de/isochrone>" params = { "locations": [{"lon": lng, "lat": lat}], "contours": [{"time": i} for i in time_steps], "costing": costing, "polygons": 1, } response = requests.post(url, json=params) if response: result = response.json() if 'error_code' not in result.keys(): df = gpd.GeoDataFrame.from_features(result) df['name'] = name df['id'] = location_id return df[['name','id','geometry']] The function takes as its input a latitude and longitude value, a costing paratemeter, a location name, and a location id. The output is a dataframe which contains a Shapely polygon representing the isochrone, along with the a name and id of the location the isochrone corresponds to. We have separate columns for a location id and a location name so that we can use the id column to examine isochrones for individual restaurants and we can use the name column to look at isochrones for each of the chains. The costing parameter can take on several different values (see the API reference here), and it can be used to create “drivetimes” assuming the user is either walking, driving, or taking public transport. We create a geoDataFrame of all of the 5-minute drivetimes for our taco restaurant locations drivetimes_5_min = pd.concat([get_isochrone(row.lat, row.long, 'auto', [5], row.chain, row.id) for row in stores.select('id','chain','lat','long').collect()]) and then save it to our S3 storage for later use: drivetimes_5_min.to_csv('s3://path/drivetimes_5_min_torchys_velvet.csv', index = False) Because we are using a free API and we have to create quite a few of these isochrones, we highly recommend saving the file for later analysis. For the purposes of this blog, we have provided a ready-made isochrone file here, which we can load into Wherobots with the following snippet: sedona.read.option('header','true').format('csv') .\\ load('s3://path/drivetimes_5_min_torchys_velvet.csv') .\\ createOrReplaceTempView('drivetimes_5_min') We can now visualize our drivetime polygons in SedonaKepler. As before, we first create the map with the snippet below. map_isochrones = sedona.read.option('header','true').format('csv'). \\ load('s3://path/drivetimes_5_min_torchys_velvet.csv') isochrone_map = SedonaKepler.create_map(map_isochrones, "Isochrones", config = isochrone_map_cfg) Now, we can display the result by calling isochrone_map . The Analysis At this point, we have a collection of the Torchy’s and Velvet Taco locations in Texas, and we know the areas which are within a 5-minute drive of each location. What we want to do now is to estimate the number of potential customers that live near each of these locations, and the extent to which these populations overlap. A First Look Before we look at how these two chains might compete with each other, let’s also take a look at the extent to which restaurants within each chain might be cannibalizing each others’ sales. A quick way to do this is by using the filtering feature in Kepler to look at isochrones for a single chain: We see that locations for each chain are fairly spread out and (at least at the 5-minute drivetime level), there is not a high degree of cannibalization within each chain. Looking at the isochrones for both chains, however, we notice that Velvet Taco locations often tend to be near Torchy’s Tacos locations (or vice-versa). At this point, all we have are qualitative statements based on these maps. Next, we will show how to use H3 and existing open-source datasets to make these statements more quantitative. Estimating Cannibalization Potential As we can see by looking at the map of isochrones above, they are highly irregular polygons which have a considerable amount of overlap. In general, these polygons are not described in a ‘nice’ way by any administrative boundaries such as census block groups, census tracts, etc. Therefore, we will have to be a little creative in order to estimate the population inside them. One way of doing this using the tools provided by Apache Sedona and Wherobots is to convert these polygons to H3 hexes. We can do this with the following snippet: sedona.sql(""" SELECT ST_H3CellIds(ST_GeomFromWKT(geometry), 8, false) AS h3, name, id FROM drivetimes_5_min """).select(explode('h3'), 'name','id').withColumnRenamed('col','h3') .\\ createOrReplaceTempView('h3_isochrones') This turns our table of drivetime polygons into a table where each row represents a hexagon with sides roughly 400m long, which is a part of a drivetime polygon. We also record the chain that these hexagons are associated to (the chain that the polygon they came from belongs to). We store each hexagon in its own row because this will simplify the process of estimating population later on. Although the question of estimating population inside individual H3 hexes is also a difficult one (we will release a notebook on this soon), open-source datasets with this information are available online, and we will use one such dataset, provided by Kontur: kontur = sedona.read.option('header','true') .\\ load('s3://path/us_h3_8_pop.geojson', format="json") .\\ drop('_corrupt_record').dropna() .\\ selectExpr('CAST(CONV(properties.h3, 16, 10) AS BIGINT) AS h3', 'properties.population as population') kontur.createOrReplaceTempView('kontur') We can now enhance our h3_isochrones table with population counts for each H3 hex: sedona.sql(""" SELECT ST_H3CellIds(ST_GeomFromWKT(geometry), 8, false) AS h3, name, id FROM drivetimes_5_min """).select(explode('h3'), 'name','id').withColumnRenamed('col','h3') .\\ join(kontur, 'h3', 'left').distinct().createOrReplaceTempView('h3_isochrones') At this stage, we can also quickly compute the cannibalization potential within each chain. Using the following code, for example, we can estimate the number of people who live within a 5 minute drive of more than one Torcy’s Tacos: sedona.sql(""" SELECT ST_H3CellIds(ST_GeomFromWKT(geometry), 8, false) AS h3, name, id FROM drivetimes_5_min """).select(explode('h3'), 'name','id').withColumnRenamed('col','h3') .\\ join(kontur, 'h3', 'left').filter('name LIKE "%Torchy%"').select('h3','population') .\\ groupBy('h3').count().filter('count >= 2').join(kontur, 'h3', 'left').distinct() .\\ agg(F.sum('population')).collect()[0][0] 97903.0 We can easily change this code to compute the same information for Velvet Taco by changing filter('name LIKE "%Torchy%"') in line 4 of the above snippet to filter('name LIKE "%Velvet%"') . If we do this, we will see that 100298 people live within a 5 minute drive of more than one Velvet Taco. Thus, we see that the Torchy’s Tacos brand appears to be slightly better at avoiding canibalization among its own locations (especially given that Torchy’s Tacos has more locations than Velvet Taco). Now, we can run the following query to show the number of people in Texas who live within a 5 minutes drive of a Torchy’s Tacos: sedona.sql(""" WITH distinct_h3 (h3, population) AS ( SELECT DISTINCT h3, ANY_VALUE(population) FROM h3_isochrones WHERE name LIKE "%Torchy's%" GROUP BY h3 ) SELECT SUM(population) FROM distinct_h3 """).show() The reason we select distinct H3 hexes here is because a single hex can belong to more than one isochrone (as evidenced by the SedonaKepler visualizations above). We get the following output: +---------------+ |sum(population)| +---------------+ | 1546765.0| +---------------+ So roughly 1.5 million people in Texas live within a 5-minute drive of a Torchy’s Tacos location. Looking at our previous calculations for how many people live near more than one restaurant of the same chain, we can see that Torchy’s Tacos locations near each other cannibalize about 6.3% of the potential customers who live within 5 minutes of a Torchy’s location. Running a similar query for Velvet Taco tells us that roughly half as many people live within a 5-minute drive of a Velvet Taco: sedona.sql(""" WITH distinct_h3 (h3, population) AS ( SELECT DISTINCT h3, ANY_VALUE(population) FROM h3_isochrones WHERE name LIKE '%Velvet Taco%' GROUP BY h3 ) SELECT SUM(population) FROM distinct_h3 """).show() +---------------+ |sum(population)| +---------------+ | 750360.0| +---------------+ As before, we can also see that Velvet Taco locations near each other cannibalize about 13.4% of the potential customers who live within 5 minutes of a Velvet Taco location. Now, we can estimate the potential for cannibalization between these two chains: sedona.sql(""" WITH overlap_h3 (h3, population) AS ( SELECT DISTINCT a.h3, ANY_VALUE(a.population) FROM h3_isochrones a LEFT JOIN h3_isochrones b ON a.h3 = b.h3 WHERE a.name != b.name GROUP BY a.h3 ) SELECT sum(population) FROM overlap_h3 """).show() which gives: +---------------+ |sum(population)| +---------------+ | 415033.0| +---------------+ We can see that more than half of the people who live near a Velvet Taco location also live near a Torchy’s Tacos location and we can visualize this population overlap: isochrones_h3_map_data = sedona.sql(""" SELECT ST_H3CellIds(ST_GeomFromWKT(geometry), 8, false) AS h3, name, id FROM drivetimes_5_min """).select(explode('h3'), 'name','id').withColumnRenamed('col','h3') .\ join(kontur, 'h3', 'left').select('name','population',array('h3')).withColumnRenamed('array(h3)','h3').selectExpr('name','population','ST_H3ToGeom(h3)[0] AS geometry') isochrones_h3_map = SedonaKepler.create_map(isochrones_h3_map_data, 'Isochrones in H3', config = isochrones_h3_map_cfg) Create a Free Account Get Started Key takeawaysThe notebook studies retail cannibalization between Torchy’s Tacos and Velvet Taco in Texas using Overture Places, 5-minute drive isochrones from the Valhalla API, H3, and Kontur population.Drive-time polygons are converted to H3 cells at resolution 8 (hex sides roughly 400 meters) and joined to Kontur US H3-8 population.About 1,546,765 people live within a 5-minute drive of a Torchy’s; about 750,360 live within 5 minutes of a Velvet Taco.97,903 people live within 5 minutes of more than one Torchy’s (about 6.3% of Torchy’s 5-minute population). 100,298 people are within 5 minutes of more than one Velvet Taco (about 13.4%).415,033 people live in the overlap of both chains’ 5-minute isochrones—more than half of Velvet Taco’s 5-minute population also lives near a Torchy’s.
Processing A Billion Aircraft Observations And Combining With Weather Data Using Apache Sedona On Wherobots Cloud Posted on May 16, 2024October 3, 2026 by Ben Pruden In this post we take a look at combining aircraft flight data with weather data to identify potentially impacted flights using public ADS-B data and weather data from NOAA. Along the way we introduce the types of spatial queries available with WherobotsDB including spatial range queries, spatial kNN, and spatial joins using Spatial SQL, geospatial visualization options, working with various file formats including GeoParquet, GeoTiff, Shapefiles, and more. We’ll be running this example in Wherobots Cloud, where you can create a free account to follow along. See the Getting Started With Wherobots Cloud series for an overview of Wherobots Cloud. Flight Trace Data The ADS-B flight observation data we’ll be using comes from ADSB.lol which provides publicly accessible daily updates of crowd-sourced flight tracking data. I adapted some great examples for downloading, parsing and enriching this dataset from Mark Litwintschik’s blog post and extracted 1 month of data which resulted in just over 1.1 billion observations. I extended Mark’s approach a bit and saved this data as GeoParquet, partitioned by S2 id. This will enable more efficient spatial queries on the dataset. First, we’ll import dependencies and configure WherobotsDB to access an S3 bucket with our example data. from sedona.spark import * config = SedonaContext.builder().config("spark.hadoop.fs.s3a.bucket.wherobots-examples.aws.credentials.provider","org.apache.hadoop.fs.s3a.AnonymousAWSCredentialsProvider").getOrCreate() sedona = SedonaContext.create(config) We load one month of world-wide aircraft trace data which is just over 1.1 billion point observations. # Load flight traces via GeoParquet traces_df = sedona.read.format("geoparquet").load("s3a://wherobots-examples/data/examples/flights/2024_s2.parquet") traces_df.createOrReplaceTempView("traces") traces_df.count() ------------------- 1103371367 Exploratory Data Analysis I followed Mark’s approach of extracting individual trace observations into their own row and we store detailed information about the trace such as altitude, squawk code, heading, etc in a struct. You can find his code in this post. Let’s view the schema of our Spatial DataFrame that we loaded via GeoParquet. traces_df.printSchema() root |-- dbFlags: long (nullable = true) |-- desc: string (nullable = true) |-- icao: string (nullable = true) |-- ownOp: string (nullable = true) |-- r: string (nullable = true) |-- reg_details: struct (nullable = true) | |-- description: string (nullable = true) | |-- iso2: string (nullable = true) | |-- iso3: string (nullable = true) | |-- nation: string (nullable = true) |-- t: string (nullable = true) |-- timestamp: string (nullable = true) |-- trace: struct (nullable = true) | |-- aircraft: struct (nullable = true) | | |-- alert: long (nullable = true) | | |-- alt_geom: long (nullable = true) | | |-- baro_rate: long (nullable = true) | | |-- category: string (nullable = true) | | |-- emergency: string (nullable = true) | | |-- flight: string (nullable = true) | | |-- geom_rate: long (nullable = true) | | |-- gva: long (nullable = true) | | |-- ias: long (nullable = true) | | |-- mach: double (nullable = true) | | |-- mag_heading: double (nullable = true) | | |-- nac_p: long (nullable = true) | | |-- nac_v: long (nullable = true) | | |-- nav_altitude_fms: long (nullable = true) | | |-- nav_altitude_mcp: long (nullable = true) | | |-- nav_heading: double (nullable = true) | | |-- nav_modes: array (nullable = true) | | | |-- element: string (containsNull = true) | | |-- nav_qnh: double (nullable = true) | | |-- nic: long (nullable = true) | | |-- nic_baro: long (nullable = true) | | |-- oat: long (nullable = true) | | |-- rc: long (nullable = true) | | |-- roll: double (nullable = true) | | |-- sda: long (nullable = true) | | |-- sil: long (nullable = true) | | |-- sil_type: string (nullable = true) | | |-- spi: long (nullable = true) | | |-- squawk: string (nullable = true) | | |-- tas: long (nullable = true) | | |-- tat: long (nullable = true) | | |-- track: double (nullable = true) | | |-- track_rate: double (nullable = true) | | |-- true_heading: double (nullable = true) | | |-- type: string (nullable = true) | | |-- version: long (nullable = true) | | |-- wd: long (nullable = true) | | |-- ws: long (nullable = true) | |-- altitude: long (nullable = true) | |-- flags: long (nullable = true) | |-- geometric_altitude: long (nullable = true) | |-- geometric_vertical_rate: long (nullable = true) | |-- ground_speed: double (nullable = true) | |-- h3_5: string (nullable = true) | |-- indicated_airspeed: long (nullable = true) | |-- lat: double (nullable = true) | |-- lon: double (nullable = true) | |-- roll_angle: double (nullable = true) | |-- source: string (nullable = true) | |-- timestamp: string (nullable = true) | |-- track_degrees: double (nullable = true) | |-- vertical_rate: long (nullable = true) |-- year: string (nullable = true) |-- geometry: geometry (nullable = true) |-- date: date (nullable = true) |-- s2: long (nullable = true) We can also inspect the first few rows of data. traces_df.show(5) +-------+--------------------+------+-----+------+--------------------+----+-------------------+--------------------+----+--------------------+----------+-------------------+ |dbFlags| desc| icao|ownOp| r| reg_details| t| timestamp| trace|year| geometry| date| s2| +-------+--------------------+------+-----+------+--------------------+----+-------------------+--------------------+----+--------------------+----------+-------------------+ | 0|PIPER PA-28R-180/...|c81c51| null|ZK-RTE|{general, NZ, NZL...|P28R|2024-01-20 00:00:00|{{null, null, nul...|null|POINT (174.722786...|2024-01-20|7854277750134145024| | 0| AIRBUS A-320|c81e2c| null|ZK-OJS|{general, NZ, NZL...|A320|2024-01-15 00:00:00|{{0, 36400, 0, A3...|null|POINT (174.549826...|2024-01-15|7854277750134145024| | 0|PIPER PA-28R-180/...|c81c51| null|ZK-RTE|{general, NZ, NZL...|P28R|2024-01-20 00:00:00|{{null, null, nul...|null|POINT (174.722764...|2024-01-20|7854277750134145024| | 0| AIRBUS A-320|c81e2c| null|ZK-OJS|{general, NZ, NZL...|A320|2024-01-15 00:00:00|{{null, null, nul...|null|POINT (174.557251...|2024-01-15|7854277750134145024| | 0|PIPER PA-28R-180/...|c81c51| null|ZK-RTE|{general, NZ, NZL...|P28R|2024-01-20 00:00:00|{{0, 1600, null, ...|null|POINT (174.722589...|2024-01-20|7854277750134145024| +-------+--------------------+------+-----+------+--------------------+----+-------------------+--------------------+----+--------------------+----------+-------------------+ only showing top 5 rows We can query our aircraft traces Spatial DataFrame using Spatial SQL. To start off let’s just demonstrate selecting fields and viewing tabular results. Note that the geometry column here is a spatial type. The data we loaded was GeoParquet formatted so we don’t need to parse or cast the field to a geometry type. sedona.sql(""" SELECT desc, ownOp, geometry, trace.ground_speed, trace.altitude, trace.aircraft.flight FROM traces WHERE trace.altitude IS NOT NULL AND OwnOp IS NOT NULL AND desc IS NOT NULL AND trace.aircraft.flight IS NOT NULL LIMIT 10 """).show(truncate=False) +----------------+-------------------+-----------------------------+------------+--------+--------+ |desc |ownOp |geometry |ground_speed|altitude|flight | +----------------+-------------------+-----------------------------+------------+--------+--------+ |AIRBUS A-350-900|UMB BANK NA TRUSTEE|POINT (174.792839 -37.010239)|122.7 |400 |DAL64 | |AIRBUS A-350-900|UMB BANK NA TRUSTEE|POINT (174.782449 -37.013079)|161.6 |350 |DAL64 | |AIRBUS A-350-900|UMB BANK NA TRUSTEE|POINT (174.779046 -37.01401) |168.9 |475 |DAL64 | |AIRBUS A-350-900|UMB BANK NA TRUSTEE|POINT (174.771463 -37.016058)|165.7 |825 |DAL64 | |AIRBUS A-350-900|UMB BANK NA TRUSTEE|POINT (174.765248 -37.017746)|162.5 |1175 |DAL64 | |AIRBUS A-350-900|UMB BANK NA TRUSTEE|POINT (174.757827 -37.01976) |159.4 |1550 |DAL64 | |AIRBUS A-350-900|UMB BANK NA TRUSTEE|POINT (174.748356 -37.022296)|161.3 |1950 |DAL64 | |AIRBUS A-350-900|UMB BANK NA TRUSTEE|POINT (174.738951 -37.024796)|163.2 |2325 |DAL64 | |AIRBUS A-350-900|UMB BANK NA TRUSTEE|POINT (174.725607 -37.028395)|165.7 |2725 |DAL64 | |AIRBUS A-350-900|UMB BANK NA TRUSTEE|POINT (174.713367 -37.03184) |167.3 |3125 |DAL64 | +----------------+-------------------+-----------------------------+------------+--------+--------+ Let’s see where the aircraft activity is in the ADS-B flight observation dataset. Where are we seeing the most aircraft activity? We can use H3 hexagons to aggregate the individual traces and visualize the results using the SedonaKepler visualization integration. Here we count the number of aircraft trace observations within each level 5 H3 hexagon and create a choropleth map to visualize the spatial distribution. h3_traces_df = sedona.sql(""" SELECT COUNT(*) AS num, ST_H3ToGeom(ST_H3CellIDs(geometry, 5, false)) AS geometry, ST_H3CellIDs(geometry, 5, false) AS h3 FROM traces GROUP BY h3 """) SedonaKepler.create_map(h3_traces_df, name="Distribution of aircraft traces") Since our data source relies on crowdsourced data retrieved from ground stations, the data coverage is not evenly distributed. We are notably missing observations over the oceans and other areas that lack crowd-sourced ground sensors. This isn’t an issue for this analysis because flight disruption is most likely to occur around airports. We can filter for aircraft belonging to a specific airline using the ownOp column and also filter for traces that were observed while the aircraft was part of a named flight. This can help give us some insight into the routes flown by this airline. jbu_h3_traces_df = sedona.sql(""" SELECT COUNT(*) AS num, ST_H3ToGeom(ST_H3CellIDs(geometry, 5, false)) AS geometry, ST_H3CellIDs(geometry, 5, false) AS h3 FROM traces WHERE ownOp = "JETBLUE AIRWAYS CORP" AND trace.aircraft.flight LIKE "JBU%" GROUP BY h3 """) SedonaKepler.create_map(jbu_h3_traces_df, name="Distribution of aircraft traces") Flight Traces – Spatial Queries And Performance WherobotsDB supports optimizations for various types of spatial queries including: Spatial range query Spatial kNN query, and Spatial joins Let’s demonstrate these types of queries and compare the performance with this dataset. Spatial Range Query Let’s say we wanted to find all aircraft traces within the bounds of a specific area, perhaps the Gulf of Mexico region. We can do this using the ST_Within Spatial SQL predicate. Here we apply this spatial predicate to filter our 1.1 billion observations to only those within the bounds of a polygon representing a specific area. This takes just over 4 seconds. gulf_filter = 'ST_Within(geometry, ST_GeomFromWKT("POLYGON ((-84.294662 29.840644, -88.952866 30.297018, -89.831772 28.767659, -94.050522 29.61167, -97.038803 27.839076, -97.917709 21.453069, -94.489975 18.479609, -86.843491 21.616579, -80.779037 24.926295, -84.294662 29.840644))"))' gulf_traces = traces_df.where(gulf_filter) %%time gulf_traces.count() -------------------- 7539347 Wall time: 4.37 s Spatial kNN Query Another type of spatial query that might be relevant for working with this dataset is a spatial kNN query, which will efficiently find the “k” nearest geometries where the value of k is specified by the user. For example, let’s find the 5 nearest aircraft traces to the Jackson Hole Airport. jackson_hole_traces = sedona.sql(""" SELECT desc, ownOp, trace.aircraft.flight, ST_DistanceSphere(ST_Point(-110.7376, 43.6088), geometry) AS distance FROM traces ORDER BY distance ASC LIMIT 5 """) %%time jackson_hole_traces.show() ---------------------------- Wall time: 11.1 s +--------------------+--------------------+--------+------------------+ | desc| ownOp| flight| distance| +--------------------+--------------------+--------+------------------+ |BOMBARDIER BD-700...| FLEXJET LLC| null| 76.92424103100763| | AIRBUS A-320| UNITED AIRLINES INC| null|161.98830220247405| |CESSNA 700 CITATI...|AH CAPITAL MANAGE...| null|180.05761405860687| |CESSNA 750 CITATI...| 2FLIU LLC| null|187.01797656145476| |DASSAULT FALCON 2000|NEXTERA ENERGY EQ...|N58MW |204.06240031528034| In this case our query takes 11 seconds against a 1.1 billion row dataset. Historical Impact of Severe Weather Events On Flights An important analysis related to flight routes is the impact of severe weather events. To enable this type of analysis Wherobots makes public a dataset of severe weather events as part of the Wherobots Open Data Catalog. The dataset has over 8 million weather events reported by NOAA in the US from 2016-2022. weather = sedona.table("wherobots_pro_data.weather.weather_events") weather.show() +--------+-------------+--------+-------------------+-------------------+-----------------+----------+-----------+-----------+-----------+----+--------------+-----+-------+--------------------+ | EventId| Type|Severity| StartTime(UTC)| EndTime(UTC)|Precipitation(in)| TimeZone|AirportCode|LocationLat|LocationLng|City| County|State|ZipCode| geometry| +--------+-------------+--------+-------------------+-------------------+-----------------+----------+-----------+-----------+-----------+----+--------------+-----+-------+--------------------+ |W-237263| Cold| Severe|2016-01-01 08:57:00|2016-01-04 15:32:00| 0.0|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237264| Rain| Light|2016-01-04 22:49:00|2016-01-04 23:30:00| 0.0|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237265| Cold| Severe|2016-01-05 00:35:00|2016-01-05 15:32:00| 0.0|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237266| Rain| Light|2016-01-05 15:32:00|2016-01-05 17:57:00| 0.2|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237267| Rain|Moderate|2016-01-05 17:57:00|2016-01-05 18:09:00| 0.14|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237268| Rain| Light|2016-01-05 18:09:00|2016-01-05 18:23:00| 0.07|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237269| Rain|Moderate|2016-01-05 18:23:00|2016-01-05 18:57:00| 0.35|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237270|Precipitation| UNK|2016-01-05 18:57:00|2016-01-05 19:41:00| 0.22|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237271| Cold| Severe|2016-01-06 00:57:00|2016-01-06 15:31:00| 0.0|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237272| Rain| Light|2016-01-06 16:57:00|2016-01-06 18:33:00| 0.16|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237273|Precipitation| UNK|2016-01-06 18:57:00|2016-01-06 19:10:00| 0.23|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237274| Rain|Moderate|2016-01-06 19:10:00|2016-01-06 19:57:00| 0.16|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237275| Rain| Light|2016-01-06 20:57:00|2016-01-06 21:27:00| 0.03|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237276| Rain|Moderate|2016-01-06 21:27:00|2016-01-06 22:57:00| 0.26|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237277| Cold| Severe|2016-01-07 00:57:00|2016-01-07 15:35:00| 0.0|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237278| Rain| Light|2016-01-07 18:10:00|2016-01-07 18:57:00| 0.0|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237279| Rain| Light|2016-01-07 23:57:00|2016-01-08 00:57:00| 0.0|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237280| Cold| Severe|2016-01-08 00:57:00|2016-01-08 15:29:00| 0.0|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237281| Cold| Severe|2016-01-09 00:57:00|2016-01-11 14:32:00| 0.0|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| |W-237282| Cold| Severe|2016-01-12 00:57:00|2016-01-12 15:57:00| 0.0|US/Pacific| KNSI| 33.2338| -119.4559|null|Ventura County| CA| null|POINT (-119.4559 ...| +--------+-------------+--------+-------------------+-------------------+-----------------+----------+-----------+-----------+-----------+----+--------------+-----+-------+--------------------+ only showing top 20 rows Severe Weather Events Let’s start our weather event analysis by finding areas that are commonly susceptible to severe weather events. We’ll again use H3 hexagons to aggregate the data, this time counting the number of severe weather events that occur within each hexagon. We’ll visualize the results as a choropleth to see areas where weather events frequently occur. weather_h3_df = sedona.sql(""" SELECT COUNT(*) AS num, ST_H3ToGeom(ST_H3CellIDs(ST_Buffer(geometry, 0.1), 5, false)) AS geometry, ST_H3CellIDs(ST_Buffer(geometry, 0.1), 5, false) AS h3 FROM wherobots_pro_data.weather.weather_events WHERE Severity = "Severe" GROUP BY h3 """) weather_h3_df.createOrReplaceTempView("weather_h3") SedonaKepler.create_map(weather_h3_df.cache(), name="Severe weather events") Historical Severe Weather Aircraft Trace Analysis Now that we’ve determined which areas are susceptible to frequent severe weather events, we can analyze flight trace observations in these areas to determine flights that may be susceptible to frequent severe weather events using a spatial join. Note that we’re joining a table with 1.1 billion observations against our weather event data, which takes just under 13 seconds in Wherobots Cloud. severe_weather_flight_traces_df= sedona.sql(""" SELECT traces.geometry AS geometry, desc, ownOp, trace.aircraft.flight AS flight, num AS eventCount FROM traces JOIN weather_h3 WHERE ST_Contains(weather_h3.geometry, traces.geometry) AND trace.aircraft.flight IS NOT null """) %%time severe_weather_flight_traces_df.count() ---------------------------- Wall time: 12.7 s 53788395 Using the flight number we can calculate an “intensity” metric that attempts to measure the potential impact of a flight to be disrupted by a severe weather event using the historical weather data. Here we calculate this intensity metric across all traces to find the most potentially impacted flights. most_impacted_flights = sedona.sql(""" SELECT sum(eventCount) AS intensity, ST_Collect(collect_list(geometry)) AS geometry, flight FROM severe_weather_traces WHERE NOT flight = " " AND NOT flight = "00000000" GROUP BY flight ORDER BY intensity DESC LIMIT 100 """) Flight Routes So far we’ve been working with individual aircraft traces (point observations), but it’s much more useful to analyze flights as routes. We can construct the individual flight routes from these aircraft traces, for a specific airline by filtering on the ownOp and trace.aircraft.flight fields. First we collect all aircraft traces that match our filter by airline, grouping by flight number. Then we can convert this collection of point geometries into a LineString geometry that represents the path the aircraft took during this flight using the ST_MakeLine Spatial SQL function. jbu_flights_df = sedona.sql(""" WITH jbu_traces AS ( SELECT any_value(trace.aircraft.flight) AS flight, collect_list(ST_Point(trace.lon, trace.lat)) AS points, any_value(CAST(date AS string)) AS date FROM traces WHERE ownOp = 'JETBLUE AIRWAYS CORP' AND trace.aircraft.flight IS NOT null GROUP BY trace.aircraft.flight, date ) SELECT flight, ST_MakeLine(points) AS geometry, date FROM jbu_traces WHERE size(points) > 100 """) jbu_flights_df.show() +--------+--------------------+----------+ | flight| geometry| date| +--------+--------------------+----------+ |JBU100 |LINESTRING (-112....|2024-01-24| |JBU101 |LINESTRING (-80.1...|2024-01-21| |JBU1016 |LINESTRING (-77.3...|2024-01-20| |JBU1044 |LINESTRING (-80.1...|2024-01-10| |JBU1053 |LINESTRING (-73.7...|2024-01-03| |JBU1055 |LINESTRING (-71.0...|2024-01-03| |JBU1060 |LINESTRING (-77.4...|2024-01-12| |JBU1060 |LINESTRING (-77.4...|2024-01-26| |JBU1067 |LINESTRING (-73.7...|2024-01-17| |JBU1068 |LINESTRING (-80.0...|2024-01-24| |JBU1114 |LINESTRING (-112....|2024-01-17| |JBU1131 |LINESTRING (-80.2...|2024-01-22| |JBU1136 |LINESTRING (-83.3...|2024-01-15| |JBU1142 |LINESTRING (-80.0...|2024-01-27| |JBU1196 |LINESTRING (-80.1...|2024-01-03| |JBU1196 |LINESTRING (-80.1...|2024-01-27| |JBU1217 |LINESTRING (-81.6...|2024-01-15| |JBU1230 |LINESTRING (-86.8...|2024-01-22| |JBU1238 |LINESTRING (-97.6...|2024-01-28| |JBU1251 |LINESTRING (-73.7...|2024-01-10| +--------+--------------------+----------+ only showing top 20 rows Identify Potential For Weather Impacted Flights Now that we’ve extracted flight routes from the individual aircraft trace point observations we can begin to determine the potential for severe weather events to impact these flights based on our historical weather event data. Creating Uniform Flight Route Segments To analyze which flights have the most potential to be impacted based on historical weather events we first create 100km segments of each flight which will allow us to identify not only entire flight routes that are potentially impacted by weather events, but specific flight segments with the most risk of disruption. To do this we’ll make use of the ability to create a user defined function in WherobotsDB using the Shapely library. We’ll use Shapely’s line interpolate functionality to create equal length segments of each flight route. from shapely.geometry import LineString, MultiLineString from sedona.sql.types import GeometryType # Define a Shapely UDF to create a series of consistent length line segments for each flight def interpolate_lines(line: GeometryType()): try: distance_delta = 1.0 # ~100km points_count = int(line.length // distance_delta) + 1 distances = (distance_delta * i for i in range(points_count)) points = [line.interpolate(distance) for distance in distances] + [line.boundary[1]] line_strings = [] for i in range(len(points) - 1): line_strings.append(LineString([points[i], points[i+1]])) return MultiLineString(line_strings) except: return None sedona.udf.register("interpolate_lines", interpolate_lines, GeometryType()) We can now call our user defined function in a spatial SQL statement and leverage WherobotsDB’s parallel execution to run this UDF. split_df = sedona.sql(""" SELECT flight, date, explode(ST_Dump(interpolate_lines(geometry))) AS geometry FROM jbu_flights """) split_df.createOrReplaceTempView("split_flights") split_df.show() +--------+----------+--------------------+ | flight| date| geometry| +--------+----------+--------------------+ |JBU100 |2024-01-24|LINESTRING (-118....| |JBU100 |2024-01-24|LINESTRING (-117....| |JBU100 |2024-01-24|LINESTRING (-116....| |JBU100 |2024-01-24|LINESTRING (-115....| |JBU100 |2024-01-24|LINESTRING (-114....| |JBU100 |2024-01-24|LINESTRING (-113....| |JBU100 |2024-01-24|LINESTRING (-112....| |JBU100 |2024-01-24|LINESTRING (-111....| |JBU100 |2024-01-24|LINESTRING (-110....| |JBU100 |2024-01-24|LINESTRING (-109....| |JBU100 |2024-01-24|LINESTRING (-108....| |JBU100 |2024-01-24|LINESTRING (-107....| |JBU100 |2024-01-24|LINESTRING (-106....| |JBU100 |2024-01-24|LINESTRING (-105....| |JBU100 |2024-01-24|LINESTRING (-104....| |JBU100 |2024-01-24|LINESTRING (-103....| |JBU100 |2024-01-24|LINESTRING (-102....| |JBU100 |2024-01-24|LINESTRING (-102....| |JBU100 |2024-01-24|LINESTRING (-101....| |JBU100 |2024-01-24|LINESTRING (-100....| +--------+----------+--------------------+ only showing top 20 rows Joining Flight Route Segments With Weather Events Now we perform a spatial join to identify flight segments that pass through areas that historically have a high frequency of severe weather events. severe_weather_split_df = sedona.sql(""" SELECT flight, date, split_flights.geometry AS geometry, num FROM split_flights JOIN weather_h3 WHERE ST_Intersects(weather_h3.geometry, split_flights.geometry) """) severe_weather_split_df.show() This gives us a measure of which individual flight segments have the highest potential for weather impact. +--------+----------+--------------------+----+ | flight| date| geometry| num| +--------+----------+--------------------+----+ |JBU2288 |2024-01-02|LINESTRING (-69.3...|1221| |JBU2288 |2024-01-02|LINESTRING (-70.3...|1774| |JBU2288 |2024-01-02|LINESTRING (-71.2...| 42| |JBU2288 |2024-01-02|LINESTRING (-71.2...|1192| |JBU2288 |2024-01-02|LINESTRING (-71.2...|1231| |JBU2288 |2024-01-02|LINESTRING (-72.2...| 672| |JBU2288 |2024-01-02|LINESTRING (-72.2...| 959| |JBU2288 |2024-01-02|LINESTRING (-73.2...| 652| |JBU2288 |2024-01-02|LINESTRING (-73.2...| 672| |JBU1053 |2024-01-03|LINESTRING (-73.7...| 416| |JBU2424 |2024-01-07|LINESTRING (-81.2...| 217| |JBU2424 |2024-01-07|LINESTRING (-81.2...| 467| |JBU2424 |2024-01-07|LINESTRING (-81.2...| 487| |JBU2424 |2024-01-07|LINESTRING (-81.2...| 611| |JBU2424 |2024-01-07|LINESTRING (-81.3...| 885| |JBU2424 |2024-01-07|LINESTRING (-81.3...| 722| |JBU2424 |2024-01-07|LINESTRING (-81.3...|2237| |JBU2424 |2024-01-07|LINESTRING (-81.3...| 722| |JBU2424 |2024-01-07|LINESTRING (-81.3...| 618| |JBU2424 |2024-01-07|LINESTRING (-81.3...| 607| +--------+----------+--------------------+----+ only showing top 20 rows Flights With The Highest Potential For Weather Impact Now we can join the flight segments and their corresponding weather intensity impact scores with the full flight routes to visualize the flights with the highest potential for weather impact based on the historical weather events. full_flight_intensity_df = sedona.sql(""" SELECT jbu_flights.flight, jbu_flights.geometry, jbu_flights.date, intensity FROM jbu_flights JOIN jbu_intensity ON jbu_flights.flight = jbu_intensity.flight AND jbu_flights.date = jbu_intensity.date """) map_view = SedonaKepler.create_map(full_flight_intensity_df.orderBy("intensity", ascending=False).limit(10), name="Flights") SedonaKepler.add_df(map_view, df=weather_h3_df, name="Weather") map_view Doppler Weather Raster Data So far we’ve been working with historical weather event data, lets now look at a simple example using more real-time weather radar data. Specifically, we’ll use NOAA’s NEXRAD Doppler weather radar surface reflectivity product to assess the impact of high precipitation at airports. We’ll load a GeoTiff raster image of the continental US to identify weather patterns, then join this raster data with the location of US airports. First, we load the raster image from S3 using the RS_FromPath function which will create an out-db raster in WherobotsDB. nexrad_raster_url = "s3://wherobots-examples/data/noaa/nexrad/n0q_202309231625.tif" nexrad_df = sedona.sql(f"SELECT RS_FromPath('{nexrad_raster_url}') AS raster") nexrad_df.createOrReplaceTempView("nexrad") nexrad_df.printSchema() root |-- raster: raster (nullable = true) Inspecting the meta data for this raster image we can see the spatial extent, spatial resolution, the CRS, and number of bands. sedona.sql("SELECT RS_MetaData(raster) from nexrad").show(truncate=False) +---------------------------------------------------------------------------+ |rs_metadata(raster) | +---------------------------------------------------------------------------+ |[-126.0025, 50.0025, 12200.0, 5400.0, 0.005, -0.005, 0.0, 0.0, 4326.0, 1.0]| +---------------------------------------------------------------------------+ We can visualize the raster using RS_AsImage. htmlDf = sedona.sql("SELECT RS_AsImage(raster, 1200) AS raster_image FROM nexrad") SedonaUtils.display_image(htmlDf) Next, we’ll load the point geometries of airports from a Shapefile. spatialRdd = ShapefileReader.readToGeometryRDD(sedona.sparkContext, "s3a://wherobots-examples/data/ne_50m_airports") poi_df = Adapter.toDf(spatialRdd, sedona) poi_df = poi_df poi_df.show(5) poi_df.createOrReplaceTempView("airports") +--------------------+---------+----------+-----+----------------+------+--------+--------+---------+--------------------+---------+ | geometry|scalerank|featurecla| type| name|abbrev|location|gps_code|iata_code| wikipedia|natlscale| +--------------------+---------+----------+-----+----------------+------+--------+--------+---------+--------------------+---------+ |POINT (113.935016...| 2| Airport|major| Hong Kong Int'l| HKG|terminal| VHHH| HKG|http://en.wikiped...| 150.000| |POINT (121.231370...| 2| Airport|major| Taoyuan| TPE|terminal| RCTP| TPE|http://en.wikiped...| 150.000| |POINT (4.76437693...| 2| Airport|major| Schiphol| AMS|terminal| EHAM| AMS|http://en.wikiped...| 150.000| |POINT (103.986413...| 2| Airport|major|Singapore Changi| SIN|terminal| WSSS| SIN|http://en.wikiped...| 150.000| |POINT (-0.4531566...| 2| Airport|major| London Heathrow| LHR| parking| EGLL| LHR|http://en.wikiped...| 150.000| +--------------------+---------+----------+-----+----------------+------+--------+--------+---------+--------------------+---------+ only showing top 5 rows And now we’re ready to join the raster radar data to our vector airport locations. This operation retrieves the pixel value (in this case a precipitation index value) from the raster at the location of each airport. We’ll filter for any airports where this index value is above 100, indicating significant precipitation. airport_precip_df = sedona.sql(""" SELECT name, iata_code, type, RS_Value(raster, geometry, 1) AS precip FROM airports JOIN nexrad """) airport_precip_df.where("precip > 100").show(truncate=False) +-------------------+---------+-----+------+ |name |iata_code|type |precip| +-------------------+---------+-----+------+ |Gen E L Logan Int'l|BOS |major|108.0 | |Durham Int'l |RDU |major|138.0 | +-------------------+---------+-----+------+ Based on this analysis we can see that Logan International Airport and Raleigh-Durham Airport are experiencing high precipitation events, which may potentially impact ground operations. Resources This blog post looked at combining aircraft trace observations from public ADS-B sensors with weather event and radar data to analyze the potential impact of severe weather events on flights. You can follow along with these examples by creating a free account on Wherobots Cloud and find the code for this example on GitHub. Interested in trying our Wherobots Pro features? Sign up here and receive $400 in credits for free. Try Wherobots for free Get Started Key takeawaysOne month of crowdsourced ADS-B traces from ADSB.lol was stored as S2-partitioned GeoParquet and loaded as 1,103,371,367 point observations in Wherobots Cloud.A Gulf of Mexico ST_Within range query returned 7,539,347 traces in 4.37 seconds. A kNN for the 5 nearest traces to Jackson Hole Airport took 11.1 seconds on the same 1.1 billion-row table.NOAA severe-weather events in the Wherobots catalog (over 8 million US events, 2016–2022) joined to the traces in 12.7 seconds, producing 53,788,395 traces in historically severe hexes.JetBlue traces were assembled into routes with ST_MakeLine, split into ~100 km segments via a Shapely UDF, and joined to severe-weather H3 cells to score flight-segment risk.A NEXRAD GeoTiff out-DB raster joined to airport points with RS_Value flagged Gen E L Logan Int'l (BOS) at precipitation index 108.0 and Durham Int'l (RDU) at 138.0 (threshold > 100).
Making Overture Maps Data More Efficient With GeoParquet And Apache Sedona Posted on April 16, 2024October 4, 2026 by Ben Pruden The latest release of Overture Maps data is now published in GeoParquet format, allowing for more efficient spatial operations when selecting subsets of the dataset and enabling interoperability with a growing ecosystem of data tooling that supports the cloud-native GeoParquet format. In this post we explore some of the motivations and advantages of publishing the Overture Maps data as GeoParquet as well as the process used by Overture to generate and publish a large scale geospatial dataset as GeoParquet using Apache Sedona. Understanding GeoParquet Built upon Apache Parquet, GeoParquet is an open-source file format specifically designed for storing geospatial data efficiently. Apache Parquet is a column-oriented data file format optimized for storing data efficiently in cloud data object stores like Amazon’s S3 service. However Parquet lacks native support for geospatial data. GeoParquet extends Parquet by defining a standard for storing geometry data and associated metadata, such as the spatial bounding box and coordinate reference system for each file. In a nutshell, GeoParquet adds the following information to the file metadata: Version: the GeoParquet specification version of this file Primary column: the main Geometry Type column should be considered when there are multiple Geometry Type columns stored in this file. Column metadata per Geometry Type column: Encoding of the geometry data. Currently only WKB format is supported. Geometry types: all possible geometry types occurred in this column, such as POINT, POLYGON, LINESTRING Coordinate reference system information, bounding box, geometry orientation, and so on In GeoParquet, all geometric data is required to be housed in Parquet’s Binary Type. This data must adhere to the encoding specifications detailed in the column metadata, which currently employs Well-Known Binary (WKB) encoding. The storage of geometries in GeoParquet using Parquet’s intrinsic binary format ensures that standard Parquet readers can seamlessly access GeoParquet files. This feature significantly enhances GeoParquet’s compatibility with existing data systems. Nevertheless, a dedicated GeoParquet reader, designed to interpret GeoParquet metadata, can utilize this information for further optimizations and to execute more advanced geospatial operations. Benefits of GeoParquet GeoParquet integrates all the advantages of Parquet, such as its compact data storage and rapid data retrieval capabilities. Additionally, GeoParquet introduces a standardized approach for storing and accessing geospatial data. This standardization significantly simplifies the process of data sharing and enhances interoperability across various systems and tools within the geospatial data processing ecosystem. Crucially, GeoParquet enables specialized GeoParquet readers to optimize geospatial queries, significantly enhancing the efficiency of data retrieval. This optimization is particularly vital for applications demanding real-time data processing and analysis. A significant contributor to this enhanced efficiency is the bounding box information (abbreviated as BBox) embedded in the GeoParquet file metadata. The BBox effectively serves as a spatial index, streamlining the process of data filtering. This is particularly useful when executing geospatially driven queries. When a spatial query is executed — for instance, searching for all points of interest within a certain area — the GeoParquet reader first checks the BBox information. If the query’s geographic area doesn’t intersect with the BBox of the data, the reader can immediately exclude that file from further consideration. This step vastly reduces the amount of data that needs to be retrieved and processed. As a testimony, with the help of Apache Sedona, Wherobots engineers created a GeoParquet version of the Overture 2023-07-26-alpha.0 data release which was not in GeoParquet format. This transformation led to a remarkable improvement in query performance. When conducting a spatial range query on the Building dataset using Apache Sedona, the process on the GeoParquet version took approximately 3 minutes, a significant reduction from the 1 hour and 30 minutes required for the same query on the original dataset. Further details and insights into this transformation are available in our blog post titled Harnessing Overture Maps Data: Apache Sedona’s Journey from Parquet to GeoParquet. The Role Of Apache Sedona The GIS software landscape offers a variety of tools capable of handling GeoParquet data, such as GeoPandas and QGIS. However, Apache Sedona stands apart in this field. As a renowned open-source cluster computing system, it is specifically tailored for processing large-scale spatial data. Apache Sedona offers a rich suite of spatial operations and analytics functions, all accessible through an intuitive Spatial SQL programming language. This unique blend of features establishes Apache Sedona as a highly effective system for geospatial data analysis. With the adoption of GeoParquet, Apache Sedona has expanded its capabilities, now able to create GeoParquet files. This integration not only enhances Apache Sedona’s data storage efficiency but also opens up broader prospects for advanced data analysis. This positions Apache Sedona as a distinct and powerful tool in the realm of geospatial data processing, differentiating it from other GIS tools by its ability to handle complex, large-scale spatial datasets efficiently. Further details can be found on Sedona’s website: https://sedona.apache.org/ Generating GeoParquet Files With Apache Sedona Leveraging Sedona’s programming guides, OMF crafts GeoParquet datasets by incorporating a geometry column sourced from its WKB geospatial data. The process involves utilizing Sedona functions, as illustrated in the python code snippet below: myDataFrame.withColumn("geometry", expr("ST_*")).selectExpr("ST_*") In its commitment to improving the spatial query experience for its customers, OMF has implemented a strategic pre-defined indexing method. This method organizes data effectively based on theme and type combinations. To further refine spatial filter pushdown performance, OMF capitalizes on Sedona’s ST_GeoHash function to generate dual GeoHash IDs for each geometry in the dataset. Coarse-Level GeoHashing (Precision Level 3): The first GeoHash, with a precision level of 3, creates cells approximately 156.5km x 156km in size. This level is ideal for grouping data based on geographic proximity. Using this GeoHash key, a large Sedona DataFrame is partitioned into numerous smaller GeoParquet files. Each file correlates to a unique 3-character GeoHash key, facilitating efficient data pruning at the Parquet file level. Importantly, this approach enables users to manually download individual GeoParquet files by referencing the GeoHash key in the file path, bypassing the need for a Parquet file reader. Fine-Level GeoHashing (Precision Level 8): Concurrently, OMF utilizes a finer GeoHash at level 8, resulting in smaller cells around 38.2m x 19m and then sorts data using this GeoHash within each GeoParquet file. This precision is used to sort data within each GeoParquet file. By doing so, OMF ensures that the data within each Parquet row-group or data page is organized based on spatial proximity. This organization is crucial for enabling efficient data pruning, particularly on the bbox column. The bbox column, stored in a native Parquet Struct Type, consists of four Double Type sub-columns: minx, maxx, miny, and maxy. It allows even standard Parquet readers to conduct pruning operations effectively. This capability is particularly beneficial as it complements the BBox information found in the GeoParquet file metadata. While the BBox metadata provides an overall spatial reference for the file, the bbox column offers a more granular level of spatial data, enabling granular data pruning at the Parquet row-group level. This dual GeoHash approach serves three core purposes: It allows users to download individual GeoParquet files from OMF without requiring a GeoParquet or Parquet reader, simply by checking the GeoHash key. GeoParquet readers can efficiently filter files using the BBox information in each file’s metadata, enhancing query speed. The strategy also enhances the functionality of standard Parquet readers which can effectively filter files using the statistical information of the bbox column, found in each file’s column and row-group metadata, without needing to recognize the file as a GeoParquet file. Analyze OMF GeoParquet Files with Sedona on Wherobots Cloud Let’s see some of the benefits of querying a large-scale GeoParquet dataset in action by analyzing the Overture Places theme using Apache Sedona. The easiest way to get started with Apache Sedona is by using the official Apache Sedona Docker image, or by creating a free hosted notebook in Wherobots Cloud. The following command will start a local Docker container running Apache Sedona and expose a Jupyter Lab environment on localhost port 8888: docker run -p 8888:8888 -p 8080:8080 -p 8081:8081 -p 4040:4040 \ apache/sedona:latest This Jupyter Lab environment can now be accessed via a web browser at http://localhost:8888 Next, the following Python code will load the latest Overture Maps GeoParquet release (at the time of writing) in Sedona: from sedona.spark import * config = SedonaContext. \ builder(). \ config("fs.s3a.aws.credentials.provider", "org.apache.hadoop.fs.s3a.AnonymousAWSCredentialsProvider"). \ getOrCreate() sedona = SedonaContext.create(config) df = sedona.read.format("geoparquet"). \ load("s3a://overturemaps-us-west-2/release/2024-01-17-alpha.0/theme=places/type=place") We can now apply a spatial predicate to filter for points of interest within the bounds of New York City, using the following Spatial SQL functions: ST_PolygonFromEnvelope – create a polygon geometry that will serve as the area roughly representing the boundaries of New York City from a specified bounding box ST_Within – filter rows for only points of interest that fall within the boundary of that polygon geometry spatial_filter = "ST_Within(geometry, ST_PolygonFromEnvelope(-74.25909, 40.477399, -73.700181, 40.917577))" df = df.where(spatial_filter) df.createOrReplaceTempView("places") When executing this spatial filter, Sedona can leverage the BBox GeoParquet metadata in each partitioned GeoParquet file to exclude any files that do not contain points of interest within the polygon defined in the spatial predicate. This results in less data being examined and scanned by Sedona and less data that needs to be retrieved over the network improving query performance. Also, note that the geometry column is interpreted as a geometry type without the need for explicit type casting thanks to the WKB serialization as part of the GeoParquet specification. To further demonstrate the functionality of Spatial SQL with Apache Sedona, let’s query for all stadiums within the area of New York City, then find all points of interest within a short walking distance of each stadium using the following Spatial SQL functions: ST_Buffer – create a new geometry that represents a buffer around each stadium that represents a short walking distance from the stadium ST_Intersects – filter for points of interest that lie within the buffer geometry, identifying points of interest within walking distance of each stadium stadium_places = sedona.sql(""" WITH stadiums AS (SELECT * FROM places WHERE categories.main = "stadium_arena") SELECT * FROM places, stadiums WHERE ST_Intersects(places.geometry, ST_Buffer(stadiums.geometry, 0.002)) Visualizing points of interest within walking distance from each stadium. You can see more detailed examples of analyzing the Overture Places dataset using Apache Sedona in this blog post. Conclusion The Overture Maps Foundation’s adoption of GeoParquet, in collaboration with Apache Sedona, marks a significant milestone in geospatial data management. The combination of GeoParquet’s efficient storage format and Apache Sedona’s spatial analytics capabilities brings unprecedented performance and scalability to geospatial dataset publishing and analysis. This integration opens up new avenues for researchers, developers, and geospatial enthusiasts, empowering them to explore and derive insights from vast geospatial datasets more efficiently than ever before. Interested in learning more about Apache Sedona? Download this free O’Reilly book and learn practical solutions for challenges working with various type of geospatial at scale. Access the free guide Download Now Key takeawaysOverture Maps now publishes releases as GeoParquet so geometry is stored as WKB in Parquet with file-level CRS, geometry types, and bounding-box metadata for spatial filter push-down.Wherobots converted the older Overture 2023-07-26-alpha.0 (not originally GeoParquet) and saw a Building spatial range query drop from about 1 hour 30 minutes to about 3 minutes in Apache Sedona.Overture’s Sedona pipeline writes dual GeoHash keys: precision 3 (~156.5 km × 156 km) to partition files, and precision 8 (~38.2 m × 19 m) to sort rows inside each file for row-group pruning.A native Parquet struct bbox (minx, maxx, miny, maxy) lets even non-GeoParquet readers prune row groups, while GeoParquet readers can skip whole files using metadata BBoxes.The notebook example loads the 2024-01-17-alpha.0 Places theme from s3a://overturemaps-us-west-2 and filters NYC with ST_Within, then finds POIs within ST_Buffer 0.002 of stadiums.
Raster Data Analysis With Spatial SQL And Apache Sedona Posted on March 8, 2024October 4, 2026 by Ben Pruden Using Spatial SQL for large scale vector and raster data One of the strengths of Apache Sedona and Wherobots Cloud is the ability to work with large scale vector and raster geospatial data together using Spatial SQL. In this post we’ll take a look at how to get started working with raster data in Sedona using Spatial SQL and some of the use cases for raster data analysis including vector / raster join operations, zonal statistics, and using raster map algebra. What is Raster Data Raster data refers to gridded data where each pixel has one or more values associated with it (called bands) and each pixel references a geographic location. Common raster data include aerial imagery, elevation models, precipitation maps, and population distribution datasets. This image from the National Ecological Observatory Network (NEON) does a great job of introducing raster data concepts using an example of aerial imagery. We can see the value of each pixel, the coordinates of each pixel in pixel space, and how pixel values map to the color rendered when the raster is visualized. We also see the spatial resolution of the data as measured by the height of each pixel, in this case 1 meter. This means that each pixel represents 1 square meter of the Earth’s surface. When working with raster data in Sedona we store rasters and their associated metadata as rows in a table and query them using Spatial SQL. Additionally, large rasters can be tiled so that each row in the table represents a single tile of a larger raster image. One of the tools available to us when analyzing raster data is map algebra: performing calculations using the individual pixel values. We can perform map algebra operations across multiple bands in a single raster (useful when computing indexes like NDVI) and also map algebra operations across different rasters (useful for time series analysis). There are several map algebra operations available using Spatial SQL, but in the image above we can see the RS_MapAlgebra function being used to compute the difference in average surface temperature between two years. Let’s see some examples of how to get started working with raster data in Sedona. If you’d like to follow along you can create a free Wherobots Cloud account and get started today. How to Load a Single Raster Let’s start with a simple example of loading a single GeoTiff from the NEON high resolution imagery dataset. We’ll then perform some basic map algebra operations on the image and introduce some new Spatial SQL functions along the way. Sedona supports a number of raster data formats, including Arc/Info ASCII grid, NetCDF, and GeoTiff. See the raster loaders documentation page for more. Rasters can be loaded as “out-of-database” (out-db) rasters from a remote file path or “in-database” (in-db). Out-db allows for managing large raster datasets stored on cloud storage and enables efficient handling of remote data, while the in-db format is useful when the raster data needs to managed within the database because of data integrity or access efficiency requirements. In most cases we typically want to be using the out-db raster format. We can use the RS_FromPath Spatial SQL function to create an out-db raster from a remote file path. I’ve loaded the example GeoTiff we’ll be using into a public S3 bucket, so when instantiating our SedonaContext we’ll need to configure anonymous access to this bucket. config = SedonaContext.builder(). \ config("spark.hadoop.fs.s3a.bucket.wherobots-examples.aws.credentials.provider","org.apache.hadoop.fs.s3a.AnonymousAWSCredentialsProvider"). \ getOrCreate() sedona = SedonaContext.create(config) Next, we load the GeoTiff and create an out-db raster using RS_FromPath, specifying the S3 URI for the single GeoTiff. file_url = "s3://wherobots-examples/data/examples/NEON_ortho.tif" neon_df = sedona.sql(f"SELECT RS_FromPath('{file_url}') AS raster") neon_df.createOrReplaceTempView("neon") neon_df.printSchema() We can see that our DataFrame has a single column raster which is of type raster . root |-- raster: raster (nullable = true) Once we’ve loaded a raster we can inspect the metadata using the RS_MetaData Spatial SQL function. vsedona.sql("SELECT RS_MetaData(raster) from neon").show(truncate=False) This function will return the upper left coordinates of the raster (in the raster’s coordinate system’s units), the width and height of the raster (in number of pixels), the spatial resolution of each pixel (in units of the rasters’s CRS), skew or rotation of the raster (if any), the SRID of the raster’s coordinate system, and the number of bands in the raster. +------------------------------------------------------------------------+ |rs_metadata(raster) | +------------------------------------------------------------------------+ |[470000.0, 7228000.0, 1000.0, 1000.0, 1.0, -1.0, 0.0, 0.0, 32606.0, 3.0]| +------------------------------------------------------------------------+ Based on this output we can observe a few important details about the raster: the upper left coordinate is 470000, 7228000 in the raster’s CRS the raster is 1000 pixels by 1000 pixels each pixel represents 1 square meter which means the raster image represents 1 square kilometer (1000 pixels x 1 meter) the SRID is 32606, which is UTM zone 6N the raster has 3 bands, representing the red, green, and blue color spectrums There are a number of other Spatial SQL functions that can be used to inspect similar metadata values, such as RS_NumBands sedona.sql("SELECT RS_NumBands(raster) FROM neon").show() +-------------------+ |rs_numbands(raster)| +-------------------+ | 3| +-------------------+ We can visualize the raster using the RS_AsImage function. htmlDf = sedona.sql("SELECT RS_AsImage(raster, 500) FROM neon") SedonaUtils.display_image(htmlDf) To obtain pixel values at specific coordinates we can use the RS_Value or RS_Values function. Note that for rasters with multiple bands we must specify which band we are interested in. We’ll also need to specify the SRID of our geometry (in this case a Point) to match that of the raster. Here we retrieve the value of the upper left pixel for band 1, which represents the red color spectrum. sedona.sql(""" SELECT RS_Value(raster, ST_SetSRID(ST_Point(470000, 7227880), 32606), 1) AS pixel_value FROM neon """).show(truncate=False) +-----------+ |pixel_value| +-----------+ |127.0 | +-----------+ Map Algebra Operations Performing raster math is one of the most common and powerful workflows when analyzing raster data. Sedona supports a number of band-specific raster math operations such as: RS_Add, RS_Divide, RS_Mean, RS_NormalizedDifference, etc RS_MapAlgebra RS_MapAlgebra is used to apply a map algebra script on a raster or multiple rasters, for example to calculate NDVI (to identify vegetation) or AWEI (to identify bodies of water). Let’s calculate the Normalized Difference Greenness Index (NDGI) for our raster. NDGI is a metric similar to NDVI, however it can be computed using only visible spectrum (RGB) whereas NDVI uses near infrared as well. For each pixel NDGI is calculated by dividing the difference of the green and red band values by the sum of the green and red band values, resulting in a value between -1 and 1. To interpret NDGI, values greater than zero are indicative of vegetation, while values less than zero are indicative of no-vegetation. To calculate NDGI we will use the RS_MapAlgebra function which uses a map algebra scripting language called Jiffle to express the computation to be performed. Note that the result of this operation is a new raster with a single band that represents the NDGI value for each pixel. ndgi_df = sedona.sql(""" SELECT RS_MapAlgebra(raster, 'D', 'out = ((rast[1] - rast[0]) / (rast[1] + rast[0]));') AS raster FROM neon """) ndgi_df.createOrReplaceTempView("ndgi") We can confirm the NDGI values are in the expected range by using the RS_SummaryStats function which will return the count, sum, mean, standard deviation, minimum, and maximum values of a single band in the raster. sedona.sql("SELECT RS_SummaryStats(raster) AS stats FROM ndgi").show(truncate=False) +----------------------------------------------------------------------------------------------------------------------+ |stats | +----------------------------------------------------------------------------------------------------------------------+ |[1000000.0, -10863.97160595989, -0.010863971605958938, 0.045695423898689726, -0.1493212669683258, 0.18421052631578946]| +----------------------------------------------------------------------------------------------------------------------+ We can evaluate single pixel values using the RS_Value function. First, selecting a pixel from the upper left area of the image that visually appears to be dirt or gravel we can confirm the value is less than zero. sedona.sql(""" SELECT RS_Value(raster, ST_SetSRID(ST_Point(470000, 7227880), 32606), 1) AS pixel_value FROM ndgi """).show(truncate=False) +--------------------+ |pixel_value | +--------------------+ |-0.09956709956709957| +--------------------+ Next, choosing a pixel from the lower left area of the image that appears to be tree cover we can confirm the value is greater than zero. sedona.sql(""" SELECT RS_Value(raster, ST_SetSRID(ST_Point(470700, 7228000), 32606), 1) AS pixel_value FROM ndgi """).show(truncate=False) +-------------------+ |pixel_value | +-------------------+ |0.03636363636363636| +-------------------+ We can also export our new NDGI raster as a GeoTiff for further analysis, either directly in the Wherobots Notebook environment or to S3 using the RS_AsGeoTiff function. Here we save the raster as a GeoTiff in the notebook environment. binary_raster_df = sedona.sql("SELECT RS_AsGeoTiff(raster) AS raster FROM ndgi") binary_raster_df.write.format("raster").mode("overwrite").save("ndgi.tif") Loading Multiple Rasters Previously we saw loading a single raster into a single row. In a typical workflow we work with many rasters across many rows. Let’s load a raster dataset from the WorldClim historical climate dataset of precipitation from 1970-2000. Each raster file represents the average precipitation across all observed years for each month of the year. Instead of the RS_FromPath function we used previously to load the raster this time we’ll use the binary file reader to load each raster file, then use RS_FromGeoTiff to convert each row to a raster. We do this so we also have access to the file path which we’ll use to determine the month for each raster which is included in the filename. PREC_URL = "s3://wherobots-examples/data/examples/world_clim/wc2.1_10m_prec" #/wc2.1_10m_prec_01.tif rawDf = sedona.read.format("binaryFile").load(PREC_URL + "/*.tif") rawDf.createOrReplaceTempView("rawdf") We’ll use a regular expression to match on the file name and extract the month for each raster, which is the last two digits of the filename. rasterDf = sedona.sql(""" SELECT RS_FromGeoTiff(content) AS raster, Int(regexp_extract(path, '(.*)([0-9]{2}).tif', 2)) AS month FROM rawdf """) rasterDf.createOrReplaceTempView("prec") rasterDf.printSchema() root |-- raster: raster (nullable = true) |-- month: integer (nullable = true) To determine the average monthly precipitation for a certain location we can use the RS_Value function. sedona.sql(""" SELECT month, RS_Value(raster, ST_POINT(-113.9940, 46.8721)) AS avg_prec FROM prec ORDER BY month ASC """).show(truncate=False) +-----+--------+ |month|avg_prec| +-----+--------+ |1 |43.0 | |2 |27.0 | |3 |31.0 | |4 |31.0 | |5 |50.0 | |6 |49.0 | |7 |27.0 | |8 |31.0 | |9 |31.0 | |10 |29.0 | |11 |35.0 | |12 |40.0 | +-----+--------+ What is Zonal Statistics Zonal statistics involves joining vector geometries to the raster and calculating statistical or aggregating values based on the pixel values that fall within each vector geometry. Let’s load a dataset of country boundaries and calculate the average yearly precipitation for each country. To do this we’ll make use of the RS_ZonalStats function to compute the average of all precipitation values within each country’s boundaries. Note that since the RS_ZonalStats function will return a value for each county for each month we will sum these monthly averages together to end up with the yearly average for each country. world_prec_df = sedona.sql(""" SELECT sum(RS_ZonalStats(prec.raster, countries.geometry, 1, 'avg', true)) AS yearly_avg_prec, any_value(countries.geometry), countries.name FROM prec, countries GROUP BY name ORDER BY yearly_avg_prec DESC """) world_prec_df.dropna().show() +------------------+--------------------+--------------------+ | yearly_avg_prec| any_value(geometry)| name| +------------------+--------------------+--------------------+ | 4937.5|MULTIPOLYGON (((1...| Micronesia| | 3526.0|MULTIPOLYGON (((1...| Palau| | 3378.75|MULTIPOLYGON (((-...| Samoa| | 3345.6875|MULTIPOLYGON (((1...| Brunei| | 3234.474358974358|MULTIPOLYGON (((1...| Solomon Is.| | 3111.0|MULTIPOLYGON (((-...| Saint Helena| | 3008.0|MULTIPOLYGON (((-...|Wallis and Futuna...| |2988.5508720930234|MULTIPOLYGON (((1...| Papua New Guinea| | 2934.628571428572|MULTIPOLYGON (((1...| Vanuatu| | 2881.220250521921|MULTIPOLYGON (((1...| Malaysia| |2831.0000000000005|MULTIPOLYGON (((-...| Costa Rica| |2807.0000000000005|MULTIPOLYGON (((-...| Faeroe Is.| | 2733.146385760997|MULTIPOLYGON (((1...| Indonesia| |2717.5348837209303|MULTIPOLYGON (((-...| Sierra Leone| | 2675.306306306307|MULTIPOLYGON (((-...| Panama| |2656.3900709219856|POLYGON ((-11.476...| Liberia| |2631.9744668068433|MULTIPOLYGON (((-...| Colombia| |2560.3461538461534|MULTIPOLYGON (((-...| Fiji| |2555.6666666666665|MULTIPOLYGON (((6...|São Tomé and Pr...| |2525.6388261851025|MULTIPOLYGON (((1...| Philippines| +------------------+--------------------+--------------------+ only showing top 20 rows Visualizing the results of this zonal statistics analysis we can see that the countries with the highest average annual precipitation lie mostly along the equator. Tiling Rasters If we are working with large rasters we typically want to break them up into multiple tiles to improve the efficiency of raster operations. Let’s take a look at an example using a dataset from NOAA of night time visible lights. A single raster image represents the intensity of light visible at night across the entire surface of the Earth. There are several derived data products from this dataset that infer measures such as population or economic output, but we’ll use the average annual visible light for a single year. We can use the RS_FromPath function to load this raster into Sedona. f10_1993_uri = "s3://wherobots-examples/data/examples/DMSP_OLS/F101993.v4b.global.stable_lights.avg_vis.tif" f10_1993_df = sedona.sql(f"SELECT RS_FromPath('{f10_1993_uri}') AS raster") f10_1993_df.createOrReplaceTempView("f10_1993") This gives us a single row with a single raster. Let’s use the RS_TileExplode function to break our single large raster into equally sized tiles, with one tile per row. tile_df = sedona.sql("SELECT RS_TileExplode(raster, 256, 256) AS (x, y, tile) FROM f10_1993") tile_df.show(5) +---+---+--------------------+ | x| y| tile| +---+---+--------------------+ | 0| 0|OutDbGridCoverage...| | 1| 0|OutDbGridCoverage...| | 2| 0|OutDbGridCoverage...| | 3| 0|OutDbGridCoverage...| | 4| 0|OutDbGridCoverage...| +---+---+--------------------+ So far we’ve been querying our data as Sedona Spatial DataFrames, without persisting these DataFrames. However, to persist these tables we can save our DataFrames as Havasu Tables. Havasu is a spatial table format based on Apache Iceberg used by Sedona for storing geospatial data and enables efficient geospatial operations. You can read more about Havasu in this post. By default in Wherobots Cloud the wherobots catalog is configured as a namespace for creating and storing databases that will be persisted in S3 cloud storage private to our Wherobots user. Here we save our tiled raster table to a database called test_db in the wherobots catalog, thus the full name for this table is wherobots.test_db.f10_1993 sedona.sql("DROP TABLE IF EXISTS wherobots.test_db.f10_1993") tile_df.writeTo("wherobots.test_db.f10_1993").create() We can then refer to this table during queries. Here we count the number of rows in the table, which in this case is the number of tiles. sedona.table("wherobots.test_db.f10_1993").count() ----------------------- 11154 To visualize the boundaries of each tile we can use SedonaKepler, along with the RS_Envelope function which will return the bounding box of each tile. Note that each tile also has x,y coordinates which represents its position relative to other tiles. tiledMap = SedonaKepler.create_map() SedonaKepler.add_df(tiledMap, sedona.table("wherobots.test_db.f10_1993").withColumn("tile", expr("RS_Envelope(tile)")), name="tiles") We can also visualize each individual tile using the RS_AsImage function. display(HTML(sedona.table("wherobots.test_db.f10_1993").selectExpr("RS_AsImage(tile, 200)", "RS_Envelope(tile)").limit(5).toPandas().to_html(escape=False))) To inspect the observed light intensity at a certain point we first need to find the tile that contains the point we wish to inspect. We can do this using the RS_Intersects predicate to find intersecting vector and raster geometries. sedona.sql(""" WITH matched_tile AS ( SELECT * FROM wherobots.test_db.f10_1993 WHERE RS_Intersects(tile, ST_POINT(-113.9940, 46.8721)) ) SELECT RS_Value(tile, ST_POINT(-113.9940, 46.8721)) AS pixel_value FROM matched_tile """).show(truncate=False) +-----------+ |pixel_value| +-----------+ |63.0 | +-----------+ Tiled Zonal Statistics We can also perform zonal statistics operations on a tiled raster dataset. This time let’s determine which counties in the US have the most night time visible light. First, we’ll load the boundaries of all US counties using a Shapefile from Natural Earth. counties_shapefile = "s3://wherobots-examples/data/examples/natural_earth/ne_10m_admin_2_counties" spatialRDD = ShapefileReader.readToGeometryRDD(sedona, counties_shapefile) counties_df = Adapter.toDf(spatialRDD, sedona) counties_df.createOrReplaceTempView("counties") Similar to the zonal statistics example above we’ll use the RS_ZonalStats function to calculate statistics of pixel values from our raster for all pixels that are contained within a geometry, however since our raster is tiled we’ll first need to match the tile(s) that intersect with the boundary each of county. county_light_tiled_df = sedona.sql(""" WITH matched_tile AS ( SELECT tile, geometry, FIPS FROM wherobots.test_db.f10_1993, counties WHERE RS_Intersects(tile, counties.geometry) ) SELECT sum(RS_ZonalStats(matched_tile.tile, matched_tile.geometry, 'mean')) AS mean_light, any_value(matched_tile.geometry) AS geometry, FIPS FROM matched_tile GROUP BY FIPS """) county_light_tiled_df.createOrReplaceTempView("county_light_1993") county_light_tiled_df.show(100) +--------------------+--------------------+-------+ | mean_light| geometry| FIPS| +--------------------+--------------------+-------+ | 6.563388510224063|MULTIPOLYGON (((-...|US01003| | 5.303058346553825|POLYGON ((-85.422...|US01019| | 7.112594570538514|POLYGON ((-86.413...|US01021| | 2.5223492723492806|POLYGON ((-88.091...|US01025| | 9.564617731305812|POLYGON ((-85.789...|US01031| | 13.433770014555993|POLYGON ((-88.130...|US01033| | 8.240051020408188|POLYGON ((-86.370...|US01037| | 2.301078582434537|POLYGON ((-86.191...|US01039| | 0.9495387954422182|POLYGON ((-86.499...|US01041| | 8.3112128146453|POLYGON ((-85.770...|US01045| | 2.008212672420753|POLYGON ((-86.916...|US01047| | 4.339487179487201|POLYGON ((-86.699...|US01053| Visualizing the results of this analysis we can see that night time visible light intensity is largely consistent with population centers. Further analysis of this data using time series analysis could reveal patterns of growth or resource extraction activities. SedonaKepler.create_map(county_light_tiled_df, name="Night Time Visible Lights by County") Resources Here are some relevant resources that you might find helpful: Working with raster data Wherobots documentation Raster function reference documentation Jupyter notebook code on GitHub Wherobots open data raster datasets Interested in learning more about Apache Sedona? Get this hands-on guide that we’ve partnered with O’Reilly on real world examples of how to leverage Apache Sedona, along with other technologies, to unlock the potential of geospatial analytics at planetary scale. Download the guide Access Now Key takeawaysApache Sedona and Wherobots store rasters as a native column type and query them with Spatial SQL, including out-DB rasters from S3 via RS_FromPath and in-DB rasters via RS_FromGeoTiff.A NEON orthoexample GeoTiff is 1,000 × 1,000 pixels at 1 meter resolution (1 km²), SRID 32606, 3 RGB bands. RS_MapAlgebra computes NDGI from green and red; RS_SummaryStats reports 1,000,000 pixels with mean about -0.011.WorldClim 1970–2000 monthly precipitation rasters are loaded as 12 rows; RS_Value at a sample point returns monthly averages such as 43.0 mm in January.Zonal statistics with RS_ZonalStats against country polygons produce yearly average precipitation (Micronesia 4,937.5 at the top of the shown table).A global night-lights GeoTiff tiled with RS_TileExplode(256, 256) wrote 11,154 tiles to a Havasu table (wherobots.test_db.f10_1993), then RS_Intersects + RS_ZonalStats ranked US counties by mean light.
Working With Files – Getting Started With Wherobots Cloud Part 3 Posted on February 16, 2024October 3, 2026 by Ben Pruden This is the third post in a series that will introduce Wherobots Cloud and WherobotsDB, covering how to get started with cloud-native geospatial analytics at scale. Part 1: Wherobots Cloud Overview Part 2: The Wherobots Notebook Environment Part 3: Working With Files (this post) In the previous post in this series we saw how to access and query data using Spatial SQL in Wherobots Cloud via the Wherobots Open Data Catalog. In this post we’re going to take a look at loading and working with our own data in Wherobots Cloud as well as creating and saving data as the result of our analysis, such as the end result of a data pipeline. We will cover importing files in various formats including CSV, GeoJSON, Shapefile, and GeoTIFF in WherobotsDB, working with AWS S3 cloud object storage, and creating GeoParquet files using Apache Sedona. If you’d like to follow along you can create a free Wherobots Cloud account at cloud.wherobots.com. Loading A CSV File From A Public S3 Bucket First, let’s explore loading a CSV file from a public AWS S3 bucket. In our SedonaContext object we’ll configure the anonymous S3 authentication provider for the bucket to ensure we can access this specific S3 bucket’s contents. Most access configuration will happen in the SedonaContext configuration object in the notebook, however we can also apply these settings when creating the notebook runtime by specifying additional spark configuration in the runtime configuration UI. See this page in the documentation for more examples of configuring cloud object storage access using the SedonaContext configuration object. from sedona.spark import * config = SedonaContext.builder(). \ config("spark.hadoop.fs.s3a.bucket.wherobots-examples.aws.credentials.provider", "org.apache.hadoop.fs.s3a.AnonymousAWSCredentialsProvider"). \ getOrCreate() sedona = SedonaContext.create(config) Access to private S3 buckets can be configured by either specifying access keys or for a more secure option using an IAM role trust policy. Now that we’ve configured the anonymous S3 credentials provider we can use the S3 URI to access objects within the bucket. In this case we’ll load a CSV file of bird species observations. We’ll use the ST_Point function to convert the individual longitude and latitude columns into a singe Point geometry column. S3_CSV_URL = "s3://wherobots-examples/data/examples/birdbuddy_oct23.csv" bb_df = sedona.read.format('csv'). \ option('header', 'true'). \ option('delimiter', ','). \ option('inferSchema', 'true'). \ load(S3_CSV_URL) bb_df = bb_df.selectExpr( 'ST_Point(anonymized_longitude, anonymized_latitude) AS location', 'timestamp', 'common_name', 'scientific_name') bb_df.createOrReplaceTempView('bb') bb_df.show(truncate=False) +----------------------------+-----------------------+-------------------------+----------------------+ |location |timestamp |common_name |scientific_name | +----------------------------+-----------------------+-------------------------+----------------------+ |POINT (-118.59075 34.393112)|2023-10-01 00:00:02.415|california scrub jay |aphelocoma californica| |POINT (-118.59075 34.393112)|2023-10-01 00:00:02.415|california scrub jay |aphelocoma californica| |POINT (-118.59075 34.393112)|2023-10-01 00:00:04.544|california scrub jay |aphelocoma californica| |POINT (-118.59075 34.393112)|2023-10-01 00:00:04.544|california scrub jay |aphelocoma californica| |POINT (-118.59075 34.393112)|2023-10-01 00:00:05.474|california scrub jay |aphelocoma californica| |POINT (-118.59075 34.393112)|2023-10-01 00:00:05.474|california scrub jay |aphelocoma californica| |POINT (-118.59075 34.393112)|2023-10-01 00:00:05.487|california scrub jay |aphelocoma californica| |POINT (-118.59075 34.393112)|2023-10-01 00:00:05.487|california scrub jay |aphelocoma californica| |POINT (-120.5542 43.804134) |2023-10-01 00:00:05.931|lesser goldfinch |spinus psaltria | |POINT (-120.5542 43.804134) |2023-10-01 00:00:05.931|lesser goldfinch |spinus psaltria | |POINT (-120.5542 43.804134) |2023-10-01 00:00:06.522|lesser goldfinch |spinus psaltria | |POINT (-120.5542 43.804134) |2023-10-01 00:00:06.522|lesser goldfinch |spinus psaltria | |POINT (-120.5542 43.804134) |2023-10-01 00:00:09.113|lesser goldfinch |spinus psaltria | |POINT (-120.5542 43.804134) |2023-10-01 00:00:09.113|lesser goldfinch |spinus psaltria | |POINT (-118.59075 34.393112)|2023-10-01 00:00:09.434|california scrub jay |aphelocoma californica| |POINT (-118.59075 34.393112)|2023-10-01 00:00:09.434|california scrub jay |aphelocoma californica| |POINT (-122.8521 46.864) |2023-10-01 00:00:17.488|red winged blackbird |agelaius phoeniceus | |POINT (-122.8521 46.864) |2023-10-01 00:00:17.488|red winged blackbird |agelaius phoeniceus | |POINT (-122.2438 47.8534) |2023-10-01 00:00:18.046|chestnut backed chickadee|poecile rufescens | |POINT (-122.2438 47.8534) |2023-10-01 00:00:18.046|chestnut backed chickadee|poecile rufescens | +----------------------------+-----------------------+-------------------------+----------------------+ We can visualize a sample of this data using SedonaKepler, the Kepler GL integration for Apache Sedona. SedonaKepler.create_map(df=bb_df.sample(0.001), name="Bird Species") Uploading A GeoJSON File Using The Wherobots Cloud File Browser Next, let’s see how we can upload our own data in Wherobots Cloud using the Wherobots Cloud File Browser UI. Wherobots Cloud accounts include secure file storage in private S3 buckets specific to our user or shared with other users of our organization. There are two options for uploading our own data into Wherobots Cloud: Via the file browser UI in the Wherobots Cloud web application Using the AWS CLI by generating temporary ingest credentials We’ll explore both options, first using the file browser UI. We will upload a GeoJSON file of the US Watershed Boundary Dataset. In the Wherobots Cloud web application, navigate to the “Files” tab then select the “data” directory. Here you’ll see two folders, one with the name customer-XXXXXXXXX and the other shared. The “customer” folder is private to each Wherobots Cloud user while the “shared” folder is private to the users within your Wherobots Cloud organization. We can upload files using the “Upload” button. Once the file is uploaded we can click on the clipboard icon to the right of the filename to copy the S3 URL for this file to access it in the notebook environment. We’ll save the S3 URL for this GeoJSON file as a variable to refer to later, note that this is private to my Wherobot’s user and not publicly accessible so if you’re following along you’ll have a different URL S3_URL_JSON = "s3://wbts-wbc-m97rcg45xi/qjnq6fcbf1/data/customer-hd1rff9kg390pk/getting_started/watershed_boundaries.geojson" WherobotsDB supports native readers for many file types, including GeoJSON so we’ll specify the “geojson” format to import our watersheds data into a Spatial DataFrame. This is a multiline GeoJSON file, where the features are contained in one large single object. WherobotsDB can also handle GeoJSON files with each feature in a single line. Refer to the documentation here for more information about working with GeoJSON files in WherobotsDB. watershed_df = sedona.read.format("geojson"). \ option("multiLine", "true"). \ load(S3_URL_JSON). \ selectExpr("explode(features) as features"). \ select("features.*") watershed_df.createOrReplaceTempView("watersheds") watershed_df.printSchema() root |-- geometry: geometry (nullable = true) |-- properties: struct (nullable = true) | |-- areaacres: double (nullable = true) | |-- areasqkm: double (nullable = true) | |-- globalid: string (nullable = true) | |-- huc6: string (nullable = true) | |-- loaddate: string (nullable = true) | |-- metasourceid: string (nullable = true) | |-- name: string (nullable = true) | |-- objectid: long (nullable = true) | |-- referencegnis_ids: string (nullable = true) | |-- shape_Area: double (nullable = true) | |-- shape_Length: double (nullable = true) | |-- sourcedatadesc: string (nullable = true) | |-- sourcefeatureid: string (nullable = true) | |-- sourceoriginator: string (nullable = true) | |-- states: string (nullable = true) | |-- tnmid: string (nullable = true) |-- type: string (nullable = true) The WherobotsDB GeoJSON loader will parse GeoJSON exactly as it is stored – as a single object so we’ll want to explode the features column which will give us rows in our DataFrame containing each feature’s geometry and a struct containing the feature’s associated properties. The geometries stored in GeoJSON are loaded as geometry types so we can operate on the DataFrame without explicitly creating geometries (as we did when loading the CSV file above). Here we filter for all watersheds that intersect with California and visualize them using SedonaKepler. california_df = sedona.sql(""" SELECT geometry, properties.name, properties.huc6 FROM watersheds WHERE properties.states LIKE "%CA%" """) SedonaKepler.create_map(df=california_df, name="California Watersheds") Uploading Files To Wherobots Cloud Using The AWS CLI If we have large files or many files to upload to Wherobots Cloud, instead of uploading files through the file browser web application we can use the AWS CLI directly by generating temporary ingest credentials. After clicking “Upload” select “Create Ingest Credentials” to create temporary credentials that can be used with the AWS CLI to upload data into your private Wherobots Cloud file storage. Once we generate our credentials we’ll need to configure the AWS CLI to use these credentials by adding them to a profile in the ~/.aws/credentials file (on Mac/Linux systems) or by running the aws configure command. See this page for more information on working with the AWS CLI. [default] aws_access_key_id = ASIAUZT33PSSUFF73PQA aws_secret_access_key = O4zTVJTLNURqcFJof6F21Rh7gIQOGuqAUJNREUBa aws_session_token = IQoJb3JpZ2luX2VjEOj//////////wEaCXVzLXdlc3QtMiJIMEYCIQCcn17jQory/9dbWjoq47cnxU4lENEE6S1akq1dEQx+4AIhALuFkR/XZtCOiw/AwEKtbCpj0IjDTR24MzSPxbSVoFOkKrICCMH//////////wEQABoMMzI5ODk4NDkxMDQ1IgwK3vI/VJyUkHJMltAqhgL6dhz0ikL2kpB7fCIKE52sw5NHlmG1LfuQkmxlhWHEJHnvFd1PrYEneDBTyXbt3Mxx8HQ86/k23zePAbm3mdOyVrrd7r9nA+cPYu5Jv93aGf+3brgGd/3fRMJy6y2Lwydsfuj/3u2/c8Ox7pTJKtcJYN14C8f3BNzTqtpR3bDjpTWG2+JMGFjgOx4lf9GuuXhW39tH8qOONA/y2lRiM00/j8cVOu1AZ/R5gRqL2/fCTFdxp9oBKJHXO8RZJ2u7/H67dmgdDFcw4T/ZuIvhEOtZ9TG2Vo9Vqb4jk+tP5E0ZhAnPyjWfhAdXD8at9/4i6S2WGhsywl5fnwBLYjRFRas5nnX7yzEEMLuQvq4GOpwBpyB0/qrzwPeRhHwb/K/ipspU2DMGSL0BFg6DoEpAOct/flMMmRYTaEhWV/Igexr746Hwox7ZOdN1gCLED1+iy8R/xcASNJJ9Yt194ItTvVtT4I6NTF13Oi50t0KqLURP43t1A67YwuZiWm+V7npUyiyHezBYLzTwf/rRi17lgpYO6/NSUSjgFeOqLDob11HGywzW6ifik+y269rq Let’s upload a Shapefile of US Forest Service roads to our Wherobots Cloud file storage and then import it using WherobotsDB. The command to upload our Shapefile data will look like this (your S3 URL will be different). Note that we use the --recursive flag since we we want to upload the multiple files that make up a Shapefile. aws s3 cp --recursive roads_shapefile s3://wbts-wbc-m97rcg45xi/qjnq6fcbf1/data/customer-hd1rff9kg390pk/getting_started/roads_shapefile Back in our notebook, we can use WherobotsDB’s Shapefile Reader which will first generate a Spatial RDD from the Shapefile which we can then convert to a Spatial DataFrame. See this page in the documentation for more information about working with Shapefile data in Wherobots. S3_URL_SHAPEFILE = "s3://wbts-wbc-m97rcg45xi/qjnq6fcbf1/data/customer-hd1rff9kg390pk/getting_started/roads_shapefile" spatialRDD = ShapefileReader.readToGeometryRDD(sedona, S3_URL_SHAPEFILE) roads_df = Adapter.toDf(spatialRDD, sedona) roads_df.printSchema() We can visualize a sample of our forest service roads using SedonaKepler. SedonaKepler.create_map(df=roads_df.sample(0.1), name="Forest Service Roads") Working With Raster Data – Loading GeoTiff Raster Files So far we’ve been working with vector data: geometries and their associated properties. We can also work with raster data in Wherobots Cloud. Let’s load a GeoTiff as an out-db raster using the RS_FromPath Spatial SQL function. Once we’ve loaded the raster image we can use Spatial SQL to manipulate and work with the band data of our raster. ortho_url = "s3://wherobots-examples/data/examples/NEON_ortho.tif" ortho_df = sedona.sql(f"SELECT RS_FromPath('{ortho_url}') AS raster") ortho_df.createOrReplaceTempView("ortho") ortho_df.show(truncate=False) For example we can use the RS_AsImage function to visualize the GeoTiff, in this case an aerial image of a forest and road scene. htmlDf = sedona.sql("SELECT RS_AsImage(raster) FROM ortho") SedonaUtils.display_image(htmlDf) This image has three bands of data: red, green, and blue pixel values. sedona.sql("SELECT RS_NumBands(raster) FROM ortho").show() +-------------------+ |rs_numbands(raster)| +-------------------+ | 3| +-------------------+ We can use the RS_MapAlgebra function to calculate the normalized different greenness index (NDGI), a metric similar to the normalized difference vegetation index (NDVI) used for quantifying the health and density of vegetation and landcover. The RS_MapAlgebra function allows us to execute complex computations using values from one or more bands and also across multiple rasters if for example we had a sequence of images across time and were interested in change detection. ndgi_df = sedona.sql(""" SELECT RS_MapAlgebra(raster, 'D', 'out = (rast[1] - rast[0]) / (rast[1] + rast[0]);') AS ndgi FROM ortho """) Writing Files With Wherobots Cloud A common workflow with Wherobots Cloud is to load several files, perform some geospatial analysis and save the results as part of a larger data pipeline, often using GeoParquet. GeoParquet is a cloud-native file format that enables efficient data storage and retrieval of geospatial data. Let’s perform some geospatial analysis using the data we loaded above and save the results as GeoParquet to our Wherobots Cloud S3 bucket. We’ll perform a spatial join of our bird observations, joining with the boundaries of our watersheds and then with a GROUP BY calculating the number of bird observations in each watershed. birdshed_df = sedona.sql(""" SELECT COUNT(*) AS num, any_value(watersheds.geometry) AS geometry, any_value(watersheds.properties.name) AS name, any_value(watersheds.properties.huc6) AS huc6 FROM bb, watersheds WHERE ST_Contains(watersheds.geometry, bb.location) AND watersheds.properties.states LIKE "%CA%" GROUP BY watersheds.properties.huc6 ORDER BY num DESC """) birdshed_df.show() Here we can see the watersheds with the highest number of bird observations in our dataset. x+------+--------------------+--------------------+------+ | num| geometry| name| huc6| +------+--------------------+--------------------+------+ |138414|MULTIPOLYGON (((-...| San Francisco Bay|180500| | 75006|MULTIPOLYGON (((-...|Ventura-San Gabri...|180701| | 74254|MULTIPOLYGON (((-...|Laguna-San Diego ...|180703| | 48452|MULTIPOLYGON (((-...|Central Californi...|180600| | 33842|MULTIPOLYGON (((-...| Lower Sacramento|180201| | 20476|MULTIPOLYGON (((-...| Santa Ana|180702| | 17014|MULTIPOLYGON (((-...| San Joaquin|180400| | 15288|MULTIPOLYGON (((-...|Northern Californ...|180101| | 13636|MULTIPOLYGON (((-...| Truckee|160501| | 9964|MULTIPOLYGON (((-...|Southern Oregon C...|171003| | 6864|MULTIPOLYGON (((-...|Tulare-Buena Vist...|180300| | 5120|MULTIPOLYGON (((-...| Northern Mojave|180902| | 3660|MULTIPOLYGON (((-...| Salton Sea|181002| | 2362|MULTIPOLYGON (((-...| Carson|160502| | 1040|MULTIPOLYGON (((-...| Lower Colorado|150301| | 814|MULTIPOLYGON (((-...| Mono-Owens Lakes|180901| | 584|MULTIPOLYGON (((-...| Klamath|180102| | 516|MULTIPOLYGON (((-...| Upper Sacramento|180200| | 456|MULTIPOLYGON (((-...| North Lahontan|180800| | 436|MULTIPOLYGON (((-...|Central Nevada De...|160600| +------+--------------------+--------------------+------+ only showing top 20 rows We can also visualize this data as a choropleth using SedonaKepler. The GeoParquet data we’ll save will include the geometry of each watershed boundary, the count of bird observations, and the name and id of each watershed. We’ll save this GeoParquet file to our Wherobots private S3 bucket. Previously we accessed our S3 URL via the Wherobots File Browser UI, but we can also access this URI as an environment variable in the notebook environment. xUSER_S3_PATH = os.environ.get("USER_S3_PATH") Since this file isn’t very large we’ll save as a single partition, but typically we would want to save partitioned GeoParquet files partitioned on a geospatial index or administrative boundary. The WherobotsDB GeoParquet writer will add the geospatial metadata when saving as GeoParquet. See this page in the documentation for more information about creating GeoParquet files with Sedona. birdshed_df.repartition(1).write.mode("overwrite"). \ format("geoparquet"). \ save(USER_S3_PATH + "geoparquet/birdshed.parquet") That was a look at working with files in Wherobots Cloud. You can get started with large-scale geospatial analytics by creating a free account at cloud.wherobots.com. Please join the Wherobots & Apache Sedona Community site and let us know what you’re working on with Wherobots, what type of examples you’d like to see next, or if you have any feedback! Try the Pro Tier for Free Get Started Key takeawaysPart 3 of the getting-started series shows how to load your own files into WherobotsDB—CSV, GeoJSON, Shapefile, and GeoTIFF—and write GeoParquet results back to S3.Public S3 is configured on SedonaContext with AnonymousAWSCredentialsProvider; private buckets can use access keys or an IAM role trust policy. Uploads go through the Files UI or temporary AWS CLI ingest credentials.Each user gets a private customer-… prefix and an organization shared/ prefix. GeoJSON is read with format geojson and explode(features); Shapefiles go through ShapefileReader into a Spatial DataFrame.A CSV of Bird Buddy observations is converted with ST_Point, joined to US watershed polygons for California, and aggregated. San Francisco Bay (HUC6 180500) leads the shown table with 138,414 observations.Results are written as GeoParquet to USER_S3_PATH (notebook environment variable) using format geoparquet. Rasters load with RS_FromPath; RS_MapAlgebra computes NDGI from the NEON ortho GeoTiff.
Exploring Global Fishing Watch Public Data With Apache Sedona, Wherobots Cloud & GeoParquet Posted on January 18, 2024September 1, 2026 by Ben Pruden This post is a hands-on look at offshore ocean infrastructure and industrial vessel activity with Apache Sedona in Wherobots Cloud using data from Global Fishing Watch. We also see how GeoParquet can be used with this data to improve the efficiency of data retrieval and enable large-scale geospatial visualization. Global Fishing Watch Overview Earlier this month researchers affiliated with Global Fishing Watch published a study in the journal Nature that used machine learning and satellite imagery to reveal that 75 percent of the world’s industrial fishing vessels are not publicly tracked while 25 percent of transport and energy vessels are not publicly tracked. The researchers analyzed 2 million gigabytes of satellite radar and optical imagery to detect vessels and offshore infrastructure and then attempted to match these identified vessels with data from public ship tracking systems. You can find their study here along with links to all code and data used in the project. Global Fishing Watch is a global non-profit organization focused on creating and sharing knowledge about human activity in the oceans to ensure fair and sustainable usage of the world’s oceans. Their focus is on using big data solutions to analyze this activity (much of it from satellite imagery) and make it publicly available through a variety of open data and data products. You can find much of their open data freely available for download here. Exploring Global Fishing Watch Public Data With Apache Sedona In this post we will focus on working with the data published by Global Fishing Watch as a result of the study mentioned above which you can find in the dataset listing for “Paolo et al. (2024). Satellite mapping reveals extensive industrial activity at sea”. To follow along, first create a free account in Wherobots Cloud then download the Global Fishing Watch data linked above. The dataset for this project includes two CSV files: offshore_infrastructure_v20231106.csv 85MB – contains offshore infrastructure identified by satellite imagery and a machine learning algorithm, including the type of infrastructure (oil, wind, etc) and location industrial_vessels_v20231013.csv 6.09GB – contains vessels detected in satellite imagery and classified as fishing or non-fishing, the location of the vessel, and if the vessel was matched to public AIS records (a system for reporting and tracking ship movements) We can download these data files and upload to Wherobots Cloud via the file browser, as shown in this image (note that I’ve uploaded a number of other files from Global Fishing Watch projects): By uploading data to Wherobots Cloud this way it will be private to our user in Wherobots Cloud and hosted in AWS S3. We can click the clipboard icon to the right of each file to retrieve the S3 URL for each file. Offshore Infrastructure Let’s start by loading the offshore infrastructure data (offshore_infrastructure_v20231106.csv). First, we’ll load this into a DataFrame. offshore_infra_df = sedona.read.format('csv'). \ option('header','true').option('delimiter', ','). \ load(S3_URL_INFRASTRUCTURE) offshore_infra_df.show(5) We can view the first few rows to get a sense of what data is included. +------------+--------------+--------------+----------------+-----------------+ |structure_id|composite_date| label| lat| lon| +------------+--------------+--------------+----------------+-----------------+ | 115958| 2019-09-01|lake_maracaibo|9.71025137100954|-71.0518565233426| | 115958| 2018-12-01|lake_maracaibo|9.71025137100954|-71.0518565233426| | 115958| 2019-07-01|lake_maracaibo|9.71025137100954|-71.0518565233426| | 115958| 2018-06-01|lake_maracaibo|9.71025137100954|-71.0518565233426| | 115958| 2020-08-01|lake_maracaibo|9.71025137100954|-71.0518565233426| +------------+--------------+--------------+----------------+-----------------+ only showing top 5 rows This dataset (the smaller of the two) has 1.4 million observations. offshore_infra_df.count() ------------------------------------- 1441242 Next, we convert the string values of lat and lon to a Point geometry using the ST_POINT Spatial SQL function. Note that since the CSV format does not include a schema or types we must explicitly cast the types of any non-string values. offshore_infra_df = offshore_infra_df.selectExpr( 'ST_POINT(CAST(lon AS float), CAST(lat AS float)) AS location', 'structure_id', 'label', 'CAST(composite_date AS date) AS composite_date' ) offshore_infra_df.createOrReplaceTempView("infrastructure") offshore_infra_df.show(5) We also cast composite_date to a date type. +--------------------+------------+--------------+--------------+ | location|structure_id| label|composite_date| +--------------------+------------+--------------+--------------+ |POINT (-71.051856...| 115958|lake_maracaibo| 2019-09-01| |POINT (-71.051856...| 115958|lake_maracaibo| 2018-12-01| |POINT (-71.051856...| 115958|lake_maracaibo| 2019-07-01| |POINT (-71.051856...| 115958|lake_maracaibo| 2018-06-01| |POINT (-71.051856...| 115958|lake_maracaibo| 2020-08-01| +--------------------+------------+--------------+--------------+ only showing top 5 rows We can count the number of observations for each classification label using a GROUP BY operation. sedona.sql("SELECT COUNT(*) AS num, label FROM infrastructure GROUP BY label").show() We can see that “oil” is the most common classification for observed offshore infrastructure in this dataset. +------+--------------+ | num| label| +------+--------------+ |518138| oil| | 42358| possible_oil| | 8169| probable_wind| | 1082| possible_wind| | 23101| probable_oil| |288390|lake_maracaibo| |436659| wind| |123345| unknown| +------+--------------+ This dataset represents a time series of observations from 2017-2021 so individual pieces of infrastrucutre may be included mulitple times. Let’s see how many distinct structures are included in the data. sedona.sql(""" WITH distinct_ids AS (SELECT DISTINCT(structure_id) FROM infrastructure) SELECT COUNT(*) FROM distinct_ids """).show() -------------------- 39979 With just under 40,000 structures we should be able to visualize them using SedonaKepler. Let’s first combine some of the classification labels to simplify the visualization. distinct_infra_df = sedona.sql(""" WITH grouped AS ( SELECT structure_id, location, CASE WHEN label IN ("possible_oil", "probable_oil", "oil") THEN 'oil' WHEN label IN ("possible_wind", "probably_wind", "wind") THEN "wind" WHEN label IN ("unknown") THEN "unknown" ELSE 'other' END AS label FROM infrastructure) SELECT any_value(location) AS location, any_value(label) AS label, structure_id FROM grouped GROUP BY structure_id """) To create visualizations using SedonaKepler we can pass the DataFrame to the create_map method: SedonaKepler.create_map(distinct_infra_df, name="Offshore Infrastructure") Here we can see the distribution of offshore infrastructure structures observed in the study. Now that we have a sense of the scale of the offshore structures identified in the study, let’s explore the industrial vessel traffic dataset. Industrial Vessels First, we’ll load the industrial vessel CSV that we uploaded to Wherobots Cloud and count the number of rows. industrial_vessel_df = sedona.read.format('csv'). \ option('header','true').option('delimiter', ','). \ load(S3_URL_INDUSTRIAL) industrial_vessel_df.count() ----------------------------------- 23088625 We have just over 23 million observations. We can print the schema of the DataFrame to get a sense of the data included. industrial_vessel_df.printSchema() root |-- timestamp: string (nullable = true) |-- detect_id: string (nullable = true) |-- scene_id: string (nullable = true) |-- lat: string (nullable = true) |-- lon: string (nullable = true) |-- mmsi: string (nullable = true) |-- length_m: string (nullable = true) |-- matching_score: string (nullable = true) |-- fishing_score: string (nullable = true) |-- matched_category: string (nullable = true) |-- overpasses_2017_2021: string (nullable = true) Similarly to the previous file we loaded, all types are strings so we’ll need to cast numeric and date fields explicitly, as well as convert the individual latitude and longitude values to a point geometry. industrial_vessel_df = industrial_vessel_df.selectExpr( 'ST_Point(CAST(lon AS float), CAST(lat AS float)) AS location', 'CAST(timestamp AS Timestamp) AS timestamp', 'detect_id', 'scene_id', 'mmsi', 'CAST(length_m AS float) AS length_m', 'CAST(matching_score AS float) AS matching_score', 'CAST(fishing_score AS float) AS fishing_score', 'matched_category', 'CAST(overpasses_2017_2021 AS integer) AS overpasses_2017_2021') industrial_vessel_df.createOrReplaceTempView('iv') The data includes a matched_category field which indicates if the vessel identified was matched to the public vessel tracking system (AIS) and also if the vessel was classified as fishing or non-fishing. sedona.sql(""" SELECT COUNT(*) AS num, matched_category FROM iv GROUP BY matched_category """).show() Over 9 million vessel observations were not matched to the public vessel tracking system, indicating that many of these vessels may not have been not properly recording their location information. +--------+------------------+ |num| matched_category| +--------+------------------+ | 2864709| matched_fishing| | 803624| matched_unknown| | 9187828| unmatched| |10232464|matched_nonfishing| +--------+------------------+ Let’s visualize an aggregation of the dataset to get a sense of the spatial distribution of the data. To do this we’ll use H3 hexagons overlaid over the area covered by the vessel observations. We will count the number of vessel observations in each hexagon’s area, and then visualize the results. To generate the H3 hexagons we’ll make use of two Spatial SQL functions: ST_H3CellIDs – return an array of H3 cell ids that cover the given geometry at the specified resolution level ST_H3ToGeom – return the H3 hexagon geometry for a given array of H3 cell ids You can read more about Sedona’s suppport for H3 in this article. h3_industrial_df = sedona.sql(""" SELECT COUNT(*) AS num, ST_H3ToGeom(ST_H3CellIDs(location, 2, false)) as geometry, ST_H3CellIDs(location, 2, false) as h3 FROM iv GROUP BY h3 """) We can now visualize this aggregated data as a choropleth map using SedonaKepler. SedonaKepler.create_map(h3_industrial_df, name="Industrial Vessels" ) Note the distorted size of the hexagons in more extreme latitudes. This is an artifact of the Web Mercator map projection. When visualizing choropleth visualizations like this we typically instead want to choose a map projection that preserves the area of each hexagon so we can more easily interpret the spatial distribution of the choropleth. To do this we can convert our DataFrame to a GeoDataFrame, specifying an Equal Earth projection then using matplotlib to visualize the GeoDataFrame. h3_gdf = geopandas.GeoDataFrame(h3_industrial_df.toPandas(), geometry='geometry', crs="4326") h3_gdf = h3_gdf.to_crs(epsg=8857) ax = h3_gdf.plot( column="num", scheme="JenksCaspall", cmap="YlOrRd", legend=False, legend_kwds={"title": "Number of Observed Industrial Activities", "fontsize": 6, "title_fontsize": 6}, figsize=(24,18) ) ax.set_axis_off() Another option for visualization is to use the contextily package to include a basemap layer, but the downside is we must revert to the Web Mercator projection. Here we visualize not just the count of overall vessel observations, but specifically those unmatched to the public tracking system. This gives us an indication of areas where illicit fishing might be concentrated. h3_unmatched_gdf = geopandas.GeoDataFrame(h3_unmatched_df.toPandas(), geometry='geometry', crs="4326") h3_unmatched_gdf = h3_unmatched_gdf.to_crs(epsg=3857) ax = h3_unmatched_gdf.plot( column="num", scheme="JenksCaspall", cmap="YlOrRd", legend=True, legend_kwds={"title": "Number of Unmatched Industrial Activities", "fontsize": 6, "title_fontsize": 6}, figsize=(12,8) ) cx.add_basemap(ax,source=cx.providers.OpenStreetMap.Mapnik) ax.set_axis_off() Now that we’ve explored the data, let’s see how we can improve the performance and efficiency of working with this data as GeoParquet. GeoParquet GeoParquet is a standard that specifies how to store geospatial data in the Apache Parquet file format. The goals of GeoParquet are largely to extend the benefits of Parquet to the geospatial ecosystem and enable efficient workflows for working with geospatial attributes in columnar data. Parquet is a popular columnar file format that enables efficient data storage and efficient retrieval. Efficient data storage is achieved using various compression and encoding methods while efficient retrieval is achieved by chunking the data in row groups and storing metadata about each row group that allows for only selecting subsets of the data within the scope of a predicate (this is known as predicate pushdown) and by storing column values next to each other allowing for retrieval of only the columns necessary to complete the query. Additionally, Parquet files can be partitioned across many individual files, further allowing query engines to ignore Parquet files outside the range of the query. These features make Parquet a natural choice for efficiently storing data in cloud object stores like AWS S3. By specifying how to store geospatial data in the Parquet format, GeoParquet takes advantage of all the benefits of Parquet, with the added bonus of supporting interoperability within the geospatial data ecosystem. Let’s look at how we can take advantage of these two benefits of GeoParquet: efficient data storage and efficient retrieval. Efficient Data Storage With GeoParquet First, let’s compare the file sizes of the original dataset (CSV) with the same data in the GeoJSON and GeoParquet formats. We can save the data as GeoJSON using Sedona: industrial_vessel_df.repartition(1).write. \ mode("overwrite"). \ format("geojson"). \ save(S3_URL_DATA_ROOT + "globalfishingwatch_industrial_vessels.json") Note that we repartition the data to save as a single file. And now to create a single un-partitioned GeoParquet file for comparison: industrial_vessel_df.repartition(1).write. \ mode("overwrite"). \ format("geoparquet"). \ save(S3_URL_DATA_ROOT + "globalfishingwatch_industrial_vessels.json") The relative file sizes are: CSV: 6.1 GB GeoJSON: 9 GB GeoParquet: 2.4 GB The same data stored as GeoParquet takes up 60% less space than the CSV equivalent and 73% less storage space than the same data stored as GeoJSON. This translates to faster query times when less data is transferred over the network. Much of this efficiency is due to the GeoParquet specification that geometries are serialized using WKB, which are then encoded and compressed using the built-in compression functionality of Parquet. Efficient Retrieval With Partitioned GeoParquet A common technique when working with Parquet files is to partition the Parquet files based on the value of some attribute. When working with GeoParquet this can be done using an administrative boundary (such as by state or country). In cases where it doesn’t make sense to partition by administrative boundary, we can instead partition using a spatial index such as S2 or H3. By partitioning our data using a spatial index we can attain efficient retrieval by leveraging another feature of the GeoParquet specification: each GeoParquet file stores the bounding box of the geometries in the file in the column metadata. This means that at query time when executing a spatial filter query, the query engine can check the bounding box for each partitioned GeoParquet file and quickly exclude those falling outside the spatial filter bounds. This is known as spatial predicate pushdown. Let’s see it in action. First, we will partition our data using the S2 spatial index and save as partitioned GeoParquet across multiple files (one file for each S2 cell): industrial_vessel_df_pq = industrial_vessel_df_pq.withColumn("s2", expr("array_max(ST_S2CellIds(location, 2))")) industrial_vessel_df_pq.repartition("s2").write. \ mode("overwrite"). \ partitionBy("s2"). \ format("geoparquet"). \ save(S3_URL_DATA_ROOT + "globalfishingwatch_industrial_vessels_s2part.parquet") industrial_vessel_df_pq_part_s2 = sedona.read.format("geoparquet").load(S3_URL_DATA_ROOT + "globalfishingwatch_industrial_vessels_s2part.parquet") Now let’s compare the performance of these various file formats by executing a spatial filter query. Spatial Filter Performance Comparison Let’s query for the number of observations in the area around the Gulf of Mexico. I drew a simple polygon to capture this area, which can be represented in WKT format as: POLYGON ((-84.294662 29.840644, -88.952866 30.297018, -89.831772 28.767659, -94.050522 29.61167, -97.038803 27.839076, -97.917709 21.453069, -94.489975 18.479609, -86.843491 21.616579, -80.779037 24.926295, -84.294662 29.840644)) Next, we’ll define a spatial predicate using the ST_Within Spatial SQL function. gulf_filter = 'ST_Within(location, ST_GeomFromWKT("POLYGON ((-84.294662 29.840644, -88.952866 30.297018, -89.831772 28.767659, -94.050522 29.61167, -97.038803 27.839076, -97.917709 21.453069, -94.489975 18.479609, -86.843491 21.616579, -80.779037 24.926295, -84.294662 29.840644))"))' Now we’ll test the performance using the CSV version of the data: %%time industrial_vessel_df.where(gulf_filter).count() ---------------------------------------------------------- 217508 (16.2 s) Note that this time includes loading the entire CSV file, converting the values from string to geometry types, then performing the spatial filter without using an index. Next, we compare the performance of executing the same spatial filter but this time using our single Parquet file. industrial_vessel_df_pq = sedona.read. \ format('geoparquet'). \ load(S3_URL_DATA_ROOT + "globalfishingwatch_industrial_vessels.parquet") %%time industrial_vessel_df_pq.where(gulf_filter).count() --------------------------------------------------------------- 217508 (4.91 s) The significant improvement here was due to loading less data over the network and more efficient deserialization of the geometry values from WKB. Next, let’s compare the performance of using the partitioned GeoParquet. Sedona will be able to exclude most GeoParquet files based on their bounding box metadata significantly reducing both the data loaded over the network but also the number of rows to filter using the spatial predicate. industrial_vessel_df_pq_part_s2 = sedona.read. \ format("geoparquet"). \ load(S3_URL_DATA_ROOT + "globalfishingwatch_industrial_vessels_s2part.parquet") %%time industrial_vessel_df_pq_part_s2.where(gulf_filter).count() ------------------------------------------------------------------------ 217508 (1.87s) To summarize, the time it takes to load the data and execute the spatial filter for each file type: CSV: 16.2s Single GeoParquet file: 4.91s Partitioned GeoParquet: 1.87s Thanks to spatial predicate pushdown enabled by GeoParquet and Apache Sedona , the partitioned GeoParquet query takes 88% less time. Another exciting benefit of GeoParquet is the types of integrations throughout the geospatial data ecosystem that it enables, such as with visualization tools. Scalable Geospatial Visualization With Lonboard The GeoParquet specification is being developed alongside the GeoArrow specification, which enables efficient in-memory columnar representation of geometries and enables efficiently transporting in-memory geospatial data between tools. One thing this type of zero-copy in-process geospatial data transport can enable is efficient geospatial data visualization by efficiently transporting data to the GPU for fast visualization rendering. A current proposal to the GeoParquet specification includes GeoArrow as a serialization format option (in addition to WKB mentioned earlier). Lonboard is a relatively new Python package from the folks at Development Seed that takes advantage of GeoParquet and GeoArrow to render large scale visualizations containing millions of geometries in seconds. Since GeoArrow is not compressed it doesn’t need to be parsed before being used and the geometries can be efficiently moved to the GPU for rendering. Note that since for visualization data is moved to the web browser and the GPU of our local machine, I downloaded the GeoParquet file locally and ran this example in a local Python notebook. If running in a cloud notebook environment the data will first need to be fetched over the network so the performance impacts won’t be as significant. Let’s see how we can use Lonboard to visualize the entirety of our industrial vessel dataset (20+ million points). First, we’ll create a GeoDataFrame loading our GeoParquet data. url = "../data/industrial_vessels.parquet" gdf = gpd.read_parquet(url) We’ll visualize matched and unmatched vessels distinctly, so we’ll split them into two different GeoDataFrames and two layers in Lonboard. We’ll color matched vessels in green and unmatched in red. unmatched_gdf = gdf.loc[gdf['matched_category']=='unmatched'] matched_gdf = gdf.loc[gdf['matched_category']=='matched_fishing'] unmatched_layer = ScatterplotLayer.from_geopandas(unmatched_gdf, get_fill_color = [200, 0, 0, 200]) matched_layer = ScatterplotLayer.from_geopandas(matched_gdf, get_fill_color = [0, 200, 0, 200]) map_ = Map(layers=[unmatched_layer, matched_layer]) map_ This visualization contains millions of point and renders in seconds in our Jupyter notebook while offering smooth panning and scrolling (since all points are rendered immediately). If you’ve ever tried to render visualizations of this scale using other tools in Jupyter this is pretty impressive! Resources Video tutorial: https://youtu.be/O7h3X3cjHww?si=ilYUDrdOxvtenOmM Code: https://github.com/johnymontana/global-fishing-watch-sedona/tree/main Create a free Wherobots Cloud account: https://wherobots.services Global Fishing Watch public data: https://globalfishingwatch.org/datasets-and-code/ GeoParquet: https://geoparquet.org/ Get Started with No Commitment Try Now Key takeawaysThe post analyzes Global Fishing Watch / Paolo et al. (2024) Nature data: researchers used 2 million gigabytes of satellite radar and optical imagery and report 75% of industrial fishing vessels and 25% of transport and energy vessels are not publicly tracked.Offshore infrastructure CSV (85 MB) has 1,441,242 observations and 39,979 distinct structures; oil is the most common label (518,138 rows), then wind (436,659).Industrial vessels CSV (6.09 GB) has 23,088,625 detections: 9,187,828 unmatched to AIS, 10,232,464 matched non-fishing, 2,864,709 matched fishing, 803,624 matched unknown (2017–2021).Stored as GeoParquet the vessel table is 2.4 GB versus 6.1 GB CSV and 9 GB GeoJSON (about 60% smaller than CSV, 73% smaller than GeoJSON).A Gulf of Mexico ST_Within count (217,508 hits) took 16.2 seconds on CSV, 4.91 seconds on a single GeoParquet file, and 1.87 seconds on S2-partitioned GeoParquet (88% less time than CSV). Lonboard is used to draw 20+ million points locally.
The Wherobots Notebook Environment – Getting Started With Wherobots Cloud Part 2 Posted on January 12, 2024October 4, 2026 by Ben Pruden The Wherobots Notebook Environment is the main development interface for developers and data scientists working with WherobotsDB in Wherobots Cloud. In this post we’ll take a look at configuring and starting a notebook environment then see how to work with WherobotsDB via Python and Spatial SQL. This is the second post in a series that will introduce Wherobots Cloud, covering how to get started with cloud-native geospatial analytics at scale. Part 1: Wherobots Cloud Overview Part 2: The Wherobots Notebook Experience (this post) Part 3: Working With Files Starting The Notebook Environment After signing in to Wherobots Cloud you’ll be prompted to configure and start a notebook runtime. Anyone can create a Wherobots Cloud account and it’s free to get started. Free tier users have access to the “Tiny” runtime with a default resource configuration suitable for development and testing. Professional tier users are able to select preconfigured runtimes with more resources or create custom runtimes. We can also add additional configuration such as AWS S3 bucket credentials or additional Spark configuration before starting the runtime and also select additional Python libraries to be installed into our environment. Many PyData and geospatial Python packages are installed by default so we can typically get started with just the default runtime configuration by clicking the “Start” button. Starting the runtime creates a Jupyter notebook environment specific to our Wherobots Cloud user. For free tier users this notebook environment will run for 2 hours then shutdown after which we’ll need to restart the runtime, while professional tier users don’t have this restriction. Example Notebooks Once our notebook runtime is available we can enter the Jupyter environment by clicking “Open Notebook”. By default we’ll land on a sample notebook which introduces some basic features of SedonaDB. If you’re not familiar with Jupyter it’s an interactive development environment organized around notebooks. Notebooks are made up of cells which can be code or also include text, images, or other interactive widgets. We can run cells individually or from the menu by selecting Run >> Run All Cells to run all cells in the notebook sequentially. This initial notebook is a simple example meant to introduce some concepts and familiarize new users with working with SedonaDB. It covers: Configuring the SedonaContext to access the Wherobots Open Data Catalog Exploring the data available in the Wherobots Open Data Catalog Using spatial SQL to query for points of interest by category Visualizing the results using SedonaKepler This initial notebook is just a starting point and you can also find a number of additional example notebooks in the notebook_example directory in Jupyter. Specifically, The sedona directory includes further examples for working with SedonaDB with Python, Spatial SQL and the Overture Maps dataset sedona-example-python.ipynb – loading data from Shapefiles, performing spatial joins, and writing as GeoParquet sedona-overture-maps.ipynb – explore the Overture Maps dataset including points of interest, administrative boundaries, and road networks The havasu directory contains examples on working with the Havasu spatial table to perform ETL and data analysis using vector and raster data havasu-iceberg-geometry-etl.ipynb– creating Havasu tables, performing spatial operations, working with spatial indexes to optimize performance havasu-iceberg-raster-etl.ipynb – working with the EuroSAT raster dataset as Havasu tables, raster operations, handling CRS transforms, and benchmarking raster geometry operations havasu-iceberg-outdb-raster-etl.ipynb – demonstrates the out-db method of working with large rasters in SedonaDB, loading a large GeoTiff and splitting into tiles, joining vector data with rasters The notebook in the sedonamaps directory shows how to make use of SedonaMaps for map matching and visualizing routes sedonamaps_example.ipynb – matching noisy GPS trajectory data to OpenStreetMap road segments and visualizing the results Note: Only 1 notebook can be run at a time. If you want to run another notebook, please shut down the kernel of the current notebook first (See instructions here). Creating Your Own Notebooks & Working With Version Control Now you can of course create your own notebooks to work with spatial data in Wherobots Cloud. Once you’ve created a notebook there are a few ways to export it, for example you can download the notebook locally but a common workflow is to check notebook changes into version control and then push changes to a system like GitHub or GitLab. Our Jupyter environment has git installed so we can check our notebooks into version control from the terminal. This is also a good way to import new notebooks. For example, let’s bring in some notebooks from this repository by opening the terminal in Jupyter and running the following command: git clone https://github.com/johnymontana/30-day-map-challenge And now these notebooks are available in our Wherobots notebook environment, so as we make changes to them we can check them in to version control and push those back to our GitHub repository. Online Resources As you’re working with Wherobots Cloud be sure to join the Wherobots Community site where you can ask questions if you get stuck and also share your projects with the community. Here are some other resources you might find useful as you explore SedonaDB and Wherobots Cloud: Wherobots Online Community – Ask questions, share your projects, explore what others are working on in the community, and connect with other members of the community Wherobots YouTube Channel – Find technical tutorials, example videos, and presentations from spatial data experts on the Wherebots YouTube Channel Wherobots Documentation – The documentation includes information about how to manage your Wherobots Cloud account, how to work with data using SedonaDB, as well as reference documentation Wherobots Blog – Keep up to date with the Wherobots and Apache Sedona community including new product announcements, technical tutorials, and highlighting spatial analytics projects This was a quick introduction to the Wherobots Notebook environment. I hope you’ll enjoy working with Wherobots Cloud and I hope to see you around the Wherobots Community Site! Try Wherobots Cloud Get Started Key takeawaysThe Wherobots notebook environment is the primary IDE for WherobotsDB: a per-user Jupyter runtime started from Wherobots Cloud after you pick a size and optional Spark/S3/library config.Free-tier users get a Tiny runtime that shuts down after 2 hours. Professional-tier users can pick larger or custom runtimes without that 2-hour cap.A sample notebook on first open covers SedonaContext, the Open Data Catalog, Spatial SQL against points of interest, and SedonaKepler maps. More examples live under notebook_example (sedona, havasu, sedonamaps).Only one notebook kernel can run at a time; shut down the current kernel before starting another.git is installed in the Jupyter terminal so you can clone repositories (example: 30-day-map-challenge) and push notebook changes to GitHub or GitLab.
