From 67b23518c1fd7f0d847fb386593bd85191b89bb5 Mon Sep 17 00:00:00 2001 From: Ed Conrad Date: Tue, 27 May 2025 17:17:19 -0600 Subject: [PATCH 1/5] DuckDB to read GeoParquet and save to feature class --- geoparquet/read_and_query_geoparquet.py | 167 ++++++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100644 geoparquet/read_and_query_geoparquet.py diff --git a/geoparquet/read_and_query_geoparquet.py b/geoparquet/read_and_query_geoparquet.py new file mode 100644 index 0000000..60981d8 --- /dev/null +++ b/geoparquet/read_and_query_geoparquet.py @@ -0,0 +1,167 @@ +"""===================================================================================================================== +Name: read_and_query_geoparquet.py +Purpose: Provides a hands-on example for the blog, "How to Read and Query GeoParquet in ArcGIS Pro with DuckDB". + +Requirements: + - DuckDB v1.1 which comes with ArcGIS Pro 3.5+ by default. + +Author: Ed Conrad +Created: 5/27/2025 +=====================================================================================================================""" + +import arcpy +import os +import duckdb +from collections import namedtuple + + +def main(): + profile = os.environ['USERPROFILE'] + pro_project = os.path.join(profile, r'Documents\ArcGIS\Projects\Data_Management\Data_Management.aprx') + aprx = arcpy.mp.ArcGISProject(pro_project) + gdb = aprx.defaultGeodatabase + arcpy.env.workspace = gdb + feature_datasets = arcpy.ListDatasets() + wgs84 = arcpy.SpatialReference(4326) + fd_name = 'Colorado' + if fd_name not in feature_datasets: + print(f"Creating a feature dataset called 'Colorado' in {gdb}") + arcpy.management.CreateFeatureDataset(gdb, out_name=fd_name, spatial_reference=wgs84) + + out_path = os.path.join(gdb, fd_name) + arcpy.env.overwriteOutput = True + + # Create a DuckDB Database instance + conn = duckdb.connect() + conn.sql('install spatial;load spatial;') # Geospatial extension that adds support for working with spatial data and functions + conn.sql('install httpfs;load httpfs;') # Adds support for reading and writing files over an HTTP(S) or S3 connection + conn.sql("set s3_region='us-west-2';") # Tells DuckDB which AWS region to use when connecting to S3 buckets + conn.sql('set enable_object_cache=true;') # Improves performance by caching S3 reads locally for reuse + + # Define a BoundingBox named tuple and populate it with Colorado's bounding box values + BoundingBox = namedtuple(typename='BoundingBox', field_names=['xmin', 'xmax', 'ymin', 'ymax']) + colorado_bbox = BoundingBox(xmin=-109, xmax=-102, ymin=37, ymax=41) + colorado_bbox = BoundingBox(xmin=-108.2, xmax=-104.5, ymin=37, ymax=41) + + # Query Colorado mountain peaks from the Overture Maps dataset using DuckDB. + # Overture Maps is an open data initiative that provides geospatial data in GeoParquet format, + # which can be queried efficiently using DuckDB’s built-in Parquet and spatial extensions. + # Dataset Schema: https://docs.overturemaps.org/schema/reference/base/land/ + # Note 1: GeoParquet stores geometries in a Well-Known Binary (WKB), which means it supports all vector + # geometry types defined by the OGC Simple Features standard. + # Note 2: we leverage DuckDB's SQL dialect to select only the columns we want, passing a boundary box to limit the + # results to Colorado, and we make use of DuckDB's read_parquet() method and the DuckDB spatial extension's + # ST_AsWKB method to read the WKB + # Overture Maps regularly updates the data; + # NOTE: GeoParquet is an immutable file type, meaning new GeoParquet files are created if there's a change + release = '2025-05-21.0' + sql = f""" + SELECT + names.primary AS name, + CAST(elevation AS INT) AS elevation_meters, + CAST(elevation * 3.28084 AS INT) AS elevation_ft, + elevation * 3.28084 > 14000 AS Is_Fourteener, + ST_AsWKB(geometry) AS wkb + FROM + read_parquet( + 's3://overturemaps-us-west-2/release/{release}/theme=base/type=land/*', + filename=true, + hive_partitioning=1 + ) + WHERE + subtype = 'physical' AND class IN ('peak','volcano') AND elevation IS NOT NULL + AND elevation * 3.28084 > 14000 + AND bbox.xmin BETWEEN {colorado_bbox.xmin} AND {colorado_bbox.xmax} + AND bbox.ymin BETWEEN {colorado_bbox.ymin} AND {colorado_bbox.ymax} + """ + + # Get DuckDB Query Relation object + duck_peaks = conn.sql(sql) + + # Insert that data read from GeoParquet into an Esri Feature Class + import arcpy + print('Creating the Colorado_Fourteeners feature class...') + peaks_fc = arcpy.management.CreateFeatureclass(out_path=out_path, + out_name='Colorado_Fourteeners', + geometry_type='POINT', + spatial_reference=wgs84).getOutput(0) + arcpy.management.AddField(in_table=peaks_fc, field_name='Name', field_type='TEXT', field_length=50) + arcpy.management.AddField(in_table=peaks_fc, field_name='Elevation_meters', field_type='SHORT') + arcpy.management.AddField(in_table=peaks_fc, field_name='Elevation_ft', field_type='SHORT') + arcpy.management.AddField(in_table=peaks_fc, field_name='Is_Fourteener', field_type='SHORT') + with arcpy.da.InsertCursor(peaks_fc, field_names=['Name', 'Elevation_meters', 'Elevation_ft', 'Is_Fourteener', + 'shape@']) as iCursor: + # Fetches a single row as a tuple + raw_row = duck_peaks.fetchone() + i = 1 + while raw_row: + if i % 5 == 0: + print(f'\t-Inserted {i} peaks...') + + # Cast the tuple to list since we need to modify it + row = list(raw_row) + + # Convert the Well-Known Binary (WKB) to an Esri Geometry object + row[-1] = arcpy.FromWKB(row[-1]) + iCursor.insertRow(row) + raw_row = duck_peaks.fetchone() + i += 1 + print(f'\tSummary: inserted {i} peaks into the Colorado_Fourteeners feature class.') + + # Get Colorado Breweries! + sql = f""" + SELECT + names.primary AS Name, + addresses[1].freeform AS Address, + addresses[1].locality AS City, + addresses[1].postcode as Zip, + phones[1] as PhoneNumber, + ST_AsWKB(geometry) AS wkb + FROM + read_parquet( + 's3://overturemaps-us-west-2/release/{release}/theme=places/type=place/*', + filename=true, + hive_partitioning=1 + ) + WHERE + categories.primary = 'brewery' + AND bbox.xmin BETWEEN {colorado_bbox.xmin} AND {colorado_bbox.xmax} + AND bbox.ymin BETWEEN {colorado_bbox.ymin} AND {colorado_bbox.ymax} + """ + + # Get DuckDB Query Relation object + duck_breweries = conn.sql(sql) + + # Insert that data read from GeoParquet into an Esri Feature Class + print('Creating the Colorado_Breweries feature class...') + breweries_fc = arcpy.management.CreateFeatureclass(out_path=out_path, + out_name='Colorado_Breweries', + geometry_type='POINT', + spatial_reference=wgs84).getOutput(0) + arcpy.management.AddField(in_table=breweries_fc, field_name='Name', field_type='TEXT', field_length=100) + arcpy.management.AddField(in_table=breweries_fc, field_name='Address', field_type='TEXT', field_length=100) + arcpy.management.AddField(in_table=breweries_fc, field_name='City', field_type='TEXT', field_length=20) + arcpy.management.AddField(in_table=breweries_fc, field_name='Zip', field_type='TEXT', field_length=10) + arcpy.management.AddField(in_table=breweries_fc, field_name='Phone_Number', field_type='TEXT', field_length=12) + + with arcpy.da.InsertCursor(breweries_fc, field_names=['Name', 'Address', 'City', 'Zip', 'Phone_Number', 'shape@']) as iCursor: + # Get a single row as a tuple + raw_row = duck_breweries.fetchone() + i = 1 + while raw_row: + if i % 50 == 0: + print(f'\t-Inserted {i} breweries...') + + # Cast the tuple to list since we need to modify it + row = list(raw_row) + + # Convert the Well-Known Binary (WKB) to an Esri Geometry object + row[-1] = arcpy.FromWKB(row[-1]) + iCursor.insertRow(row) + raw_row = duck_breweries.fetchone() + i += 1 + print(f'\tSummary: inserted {i} breweries into the Colorado_Breweries feature class.') + + +if __name__ == '__main__': + main() From 5ea188e9572a3836237c7baddc0f9acbe2e92299 Mon Sep 17 00:00:00 2001 From: Ed Conrad Date: Wed, 28 May 2025 16:19:41 -0600 Subject: [PATCH 2/5] Closest Facility Analysis with Fourteeners and Breweries --- geoparquet/read_and_query_geoparquet.py | 192 ++++++++++++++++++++++-- 1 file changed, 179 insertions(+), 13 deletions(-) diff --git a/geoparquet/read_and_query_geoparquet.py b/geoparquet/read_and_query_geoparquet.py index 60981d8..3a48261 100644 --- a/geoparquet/read_and_query_geoparquet.py +++ b/geoparquet/read_and_query_geoparquet.py @@ -9,13 +9,13 @@ Created: 5/27/2025 =====================================================================================================================""" -import arcpy import os import duckdb from collections import namedtuple def main(): + import arcpy profile = os.environ['USERPROFILE'] pro_project = os.path.join(profile, r'Documents\ArcGIS\Projects\Data_Management\Data_Management.aprx') aprx = arcpy.mp.ArcGISProject(pro_project) @@ -41,9 +41,8 @@ def main(): # Define a BoundingBox named tuple and populate it with Colorado's bounding box values BoundingBox = namedtuple(typename='BoundingBox', field_names=['xmin', 'xmax', 'ymin', 'ymax']) colorado_bbox = BoundingBox(xmin=-109, xmax=-102, ymin=37, ymax=41) - colorado_bbox = BoundingBox(xmin=-108.2, xmax=-104.5, ymin=37, ymax=41) - # Query Colorado mountain peaks from the Overture Maps dataset using DuckDB. + # region Get Colorado Fourteeners! # Overture Maps is an open data initiative that provides geospatial data in GeoParquet format, # which can be queried efficiently using DuckDB’s built-in Parquet and spatial extensions. # Dataset Schema: https://docs.overturemaps.org/schema/reference/base/land/ @@ -80,9 +79,9 @@ def main(): # Insert that data read from GeoParquet into an Esri Feature Class import arcpy - print('Creating the Colorado_Fourteeners feature class...') + print('Creating the Fourteeners feature class...') peaks_fc = arcpy.management.CreateFeatureclass(out_path=out_path, - out_name='Colorado_Fourteeners', + out_name='Fourteeners', geometry_type='POINT', spatial_reference=wgs84).getOutput(0) arcpy.management.AddField(in_table=peaks_fc, field_name='Name', field_type='TEXT', field_length=50) @@ -93,7 +92,7 @@ def main(): 'shape@']) as iCursor: # Fetches a single row as a tuple raw_row = duck_peaks.fetchone() - i = 1 + i = 0 while raw_row: if i % 5 == 0: print(f'\t-Inserted {i} peaks...') @@ -106,9 +105,10 @@ def main(): iCursor.insertRow(row) raw_row = duck_peaks.fetchone() i += 1 - print(f'\tSummary: inserted {i} peaks into the Colorado_Fourteeners feature class.') + print(f'\tSummary: {i} points exist in the Fourteeners feature class.') + # endregion - # Get Colorado Breweries! + # region Get Colorado Breweries! sql = f""" SELECT names.primary AS Name, @@ -120,7 +120,7 @@ def main(): FROM read_parquet( 's3://overturemaps-us-west-2/release/{release}/theme=places/type=place/*', - filename=true, + filename=false, hive_partitioning=1 ) WHERE @@ -133,9 +133,9 @@ def main(): duck_breweries = conn.sql(sql) # Insert that data read from GeoParquet into an Esri Feature Class - print('Creating the Colorado_Breweries feature class...') + print('Creating the Breweries feature class...') breweries_fc = arcpy.management.CreateFeatureclass(out_path=out_path, - out_name='Colorado_Breweries', + out_name='Breweries', geometry_type='POINT', spatial_reference=wgs84).getOutput(0) arcpy.management.AddField(in_table=breweries_fc, field_name='Name', field_type='TEXT', field_length=100) @@ -147,7 +147,7 @@ def main(): with arcpy.da.InsertCursor(breweries_fc, field_names=['Name', 'Address', 'City', 'Zip', 'Phone_Number', 'shape@']) as iCursor: # Get a single row as a tuple raw_row = duck_breweries.fetchone() - i = 1 + i = 0 while raw_row: if i % 50 == 0: print(f'\t-Inserted {i} breweries...') @@ -160,7 +160,173 @@ def main(): iCursor.insertRow(row) raw_row = duck_breweries.fetchone() i += 1 - print(f'\tSummary: inserted {i} breweries into the Colorado_Breweries feature class.') + print(f'\tSummary: {i} points exist in the Breweries feature class.') + # endregion + + # region Get Colorado Road Segments! + # According to the docs, the definition for subtype of 'road': + # "A road segment represents a section of any kind of road, street or path, including a dedicated path for walking or cycling, but excluding a railway." + # As of 5/28/25 the names.primary is returning segment type - which differs from the docs. + sql = f""" + SELECT + names.primary as SegmentType, + ST_AsWKB(geometry) AS wkb + FROM + read_parquet( + 's3://overturemaps-us-west-2/release/{release}/theme=transportation/type=segment/*', + filename=false, + hive_partitioning=1 + ) + WHERE + subtype = 'road' + AND bbox.xmin BETWEEN {colorado_bbox.xmin} AND {colorado_bbox.xmax} + AND bbox.ymin BETWEEN {colorado_bbox.ymin} AND {colorado_bbox.ymax} + """ + + # Get DuckDB Query Relation object + duck_segments = conn.sql(sql) + + # Insert that data read from GeoParquet into an Esri Feature Class + print('Creating the Route_Segments feature class...') + segments_fc = arcpy.management.CreateFeatureclass(out_path=out_path, + out_name='Route_Segments', + geometry_type='POLYLINE', + spatial_reference=wgs84).getOutput(0) + arcpy.management.AddField(in_table=segments_fc, field_name='Segment_Type', field_type='TEXT', field_length=20) + with arcpy.da.InsertCursor(segments_fc, field_names=['shape@', 'Segment_Type']) as iCursor: + # Get a single row as a tuple + raw_row = duck_segments.fetchone() + i = 0 + while raw_row: + if i % 25_000 == 0: + print(f'\t-Inserted {i:,} segments...') + + # Cast the tuple to list since we need to modify it + row = list(raw_row) + + # Convert the Well-Known Binary (WKB) string to an Esri Geometry object + row[0] = arcpy.FromWKB(row[-1]) + iCursor.insertRow(row) + raw_row = duck_segments.fetchone() + i += 1 + print(f'\tSummary: {i:,} polylines exist in the Route_Segments feature class.') + + # region Get Colorado Road Connectors! + sql = f""" + SELECT + ST_AsWKB(geometry) AS wkb + FROM + read_parquet( + 's3://overturemaps-us-west-2/release/{release}/theme=transportation/type=connector/*', + filename=false, + hive_partitioning=1 + ) + WHERE + bbox.xmin BETWEEN {colorado_bbox.xmin} AND {colorado_bbox.xmax} + AND bbox.ymin BETWEEN {colorado_bbox.ymin} AND {colorado_bbox.ymax} + """ + + # Get DuckDB Query Relation object + duck_connectors = conn.sql(sql) + + # Insert that data read from GeoParquet into an Esri Feature Class + print('Creating the Route_Connectors feature class...') + connectors_fc = arcpy.management.CreateFeatureclass(out_path=out_path, + out_name='Route_Connectors', + geometry_type='POINT', + spatial_reference=wgs84).getOutput(0) + + with arcpy.da.InsertCursor(connectors_fc, field_names=['shape@']) as iCursor: + # Get a single row as a tuple + raw_row = duck_connectors.fetchone() + i = 0 + while raw_row: + if i % 10_000 == 0: + print(f'\t-Inserted {i:,} connectors...') + + # Cast the tuple to list since we need to modify it + row = list(raw_row) + + # Convert the Well-Known Binary (WKB) string to an Esri Geometry object + row[0] = arcpy.FromWKB(row[-1]) + iCursor.insertRow(row) + raw_row = duck_connectors.fetchone() + i += 1 + print(f'\tSummary: {i:,} points exist in the Route_Connectors feature class.') + # endregion + + # Close DuckDB connection + conn.close() + + # region Network Analyst section + # Check out Network Analyst license if available. Fail if the Network Analyst license is not available. + if arcpy.CheckExtension('Network') == 'Available': + arcpy.CheckOutExtension('Network') + else: + raise arcpy.ExecuteError('Network Analyst Extension license is not available.') + + import arcpy.na + print('Creating Route_Network...') + arcpy.na.CreateNetworkDataset(feature_dataset=out_path, + out_name='Route_Network', + source_feature_class_names='Route_Segments;Route_Connectors', + elevation_model='NO_ELEVATION') + network = os.path.join(out_path, 'Route_Network') + + print('Creating Closest Facility Analysis Layers') + closest_facility_lyr = arcpy.na.MakeClosestFacilityAnalysisLayer(network_data_source=network, + layer_name='ClosestBreweriesToFourteeners', + travel_direction='TO_FACILITIES', + number_of_facilities_to_find=3).getOutput(0) + print('Building Network...') + arcpy.na.BuildNetwork(network) + print('Network built!') + + # Add Breweries as the "Facilities" + print('Adding Breweries as the Facilities') + arcpy.na.AddLocations( + in_network_analysis_layer=closest_facility_lyr, + sub_layer='Facilities', + in_table=breweries_fc, + search_tolerance='100 Meters', # NOTE: value determined after some trial-and-error + match_type='MATCH_TO_CLOSEST', + append='APPEND', + snap_to_position_along_network='SNAP', + snap_offset='0 Meters' + ) + + # Add Fourteeners as the "Incidents" + print('Adding Fourteeners as the Incidents') + arcpy.na.AddLocations( + in_network_analysis_layer=closest_facility_lyr, + sub_layer='Incidents', + in_table=peaks_fc, + search_tolerance='1500 Meters', # NOTE: value determined after some trial-and-error + match_type='MATCH_TO_CLOSEST', + append='APPEND', + snap_to_position_along_network='SNAP', + snap_offset='0 Meters' + ) + + print('\n\nNow you need to make a manual property change to the Network Dataset that was just created.\n\n' + 'Instructions:\n' + f'-Expand {gdb}\n' + '-Expand the "Colorado" feature dataset\n' + '-Right-click the Route_Network dataset that was just created\n' + '-Select Properties -> Source Settings -> Group Connectivity\n' + '-Change the policy for Route_Connectors from "Honor" to "Override".\n\n') + print('After making this change, uncomment the code below and run it to rebuild the Network to incorporate the change and solve.') + + # print('Rebuilding Network...') + # arcpy.na.BuildNetwork(network) + # print('Network rebuilt!') + # + # print('Solving nearest 3 breweries to each fourteener') + # arcpy.na.Solve(closest_facility_lyr) + # print('Routes from each Fourteener to nearest 3 breweries created!') + # if arcpy.CheckExtension('Network') == 'Available': + # arcpy.CheckInExtension('Network') + # endregion if __name__ == '__main__': From 462e1587ce1c2b276167d3da3cdaadeecc01c691 Mon Sep 17 00:00:00 2001 From: Ed Conrad Date: Wed, 28 May 2025 16:32:27 -0600 Subject: [PATCH 3/5] simplify read_parquet() function calls --- geoparquet/read_and_query_geoparquet.py | 26 ++++++------------------- 1 file changed, 6 insertions(+), 20 deletions(-) diff --git a/geoparquet/read_and_query_geoparquet.py b/geoparquet/read_and_query_geoparquet.py index 3a48261..b0fa9e0 100644 --- a/geoparquet/read_and_query_geoparquet.py +++ b/geoparquet/read_and_query_geoparquet.py @@ -51,6 +51,8 @@ def main(): # Note 2: we leverage DuckDB's SQL dialect to select only the columns we want, passing a boundary box to limit the # results to Colorado, and we make use of DuckDB's read_parquet() method and the DuckDB spatial extension's # ST_AsWKB method to read the WKB + # Docs for read_parquet() https://duckdb.org/docs/stable/data/parquet/overview#read_parquet-function + # Docs for ST_AsWKB() https://duckdb.org/docs/stable/core_extensions/spatial/functions#st_aswkb # Overture Maps regularly updates the data; # NOTE: GeoParquet is an immutable file type, meaning new GeoParquet files are created if there's a change release = '2025-05-21.0' @@ -62,11 +64,7 @@ def main(): elevation * 3.28084 > 14000 AS Is_Fourteener, ST_AsWKB(geometry) AS wkb FROM - read_parquet( - 's3://overturemaps-us-west-2/release/{release}/theme=base/type=land/*', - filename=true, - hive_partitioning=1 - ) + read_parquet('s3://overturemaps-us-west-2/release/{release}/theme=base/type=land/*') WHERE subtype = 'physical' AND class IN ('peak','volcano') AND elevation IS NOT NULL AND elevation * 3.28084 > 14000 @@ -118,11 +116,7 @@ def main(): phones[1] as PhoneNumber, ST_AsWKB(geometry) AS wkb FROM - read_parquet( - 's3://overturemaps-us-west-2/release/{release}/theme=places/type=place/*', - filename=false, - hive_partitioning=1 - ) + read_parquet('s3://overturemaps-us-west-2/release/{release}/theme=places/type=place/*') WHERE categories.primary = 'brewery' AND bbox.xmin BETWEEN {colorado_bbox.xmin} AND {colorado_bbox.xmax} @@ -172,11 +166,7 @@ def main(): names.primary as SegmentType, ST_AsWKB(geometry) AS wkb FROM - read_parquet( - 's3://overturemaps-us-west-2/release/{release}/theme=transportation/type=segment/*', - filename=false, - hive_partitioning=1 - ) + read_parquet('s3://overturemaps-us-west-2/release/{release}/theme=transportation/type=segment/*') WHERE subtype = 'road' AND bbox.xmin BETWEEN {colorado_bbox.xmin} AND {colorado_bbox.xmax} @@ -216,11 +206,7 @@ def main(): SELECT ST_AsWKB(geometry) AS wkb FROM - read_parquet( - 's3://overturemaps-us-west-2/release/{release}/theme=transportation/type=connector/*', - filename=false, - hive_partitioning=1 - ) + read_parquet('s3://overturemaps-us-west-2/release/{release}/theme=transportation/type=connector/*') WHERE bbox.xmin BETWEEN {colorado_bbox.xmin} AND {colorado_bbox.xmax} AND bbox.ymin BETWEEN {colorado_bbox.ymin} AND {colorado_bbox.ymax} From 3d0ada5e4cba1a5c89c0c651986fd429bb8acc4b Mon Sep 17 00:00:00 2001 From: Ed Conrad Date: Wed, 28 May 2025 18:01:46 -0600 Subject: [PATCH 4/5] added logic to print out summary of the results --- geoparquet/read_and_query_geoparquet.py | 74 +++++++++++++++++++------ 1 file changed, 56 insertions(+), 18 deletions(-) diff --git a/geoparquet/read_and_query_geoparquet.py b/geoparquet/read_and_query_geoparquet.py index b0fa9e0..9693020 100644 --- a/geoparquet/read_and_query_geoparquet.py +++ b/geoparquet/read_and_query_geoparquet.py @@ -10,9 +10,12 @@ =====================================================================================================================""" import os -import duckdb from collections import namedtuple +import arcpy.da +import arcpy.management +import duckdb + def main(): import arcpy @@ -227,7 +230,7 @@ def main(): raw_row = duck_connectors.fetchone() i = 0 while raw_row: - if i % 10_000 == 0: + if i % 25_000 == 0: print(f'\t-Inserted {i:,} connectors...') # Cast the tuple to list since we need to modify it @@ -294,26 +297,61 @@ def main(): snap_offset='0 Meters' ) - print('\n\nNow you need to make a manual property change to the Network Dataset that was just created.\n\n' - 'Instructions:\n' - f'-Expand {gdb}\n' - '-Expand the "Colorado" feature dataset\n' - '-Right-click the Route_Network dataset that was just created\n' - '-Select Properties -> Source Settings -> Group Connectivity\n' - '-Change the policy for Route_Connectors from "Honor" to "Override".\n\n') + print('\n\nNow you need to make a manual property change to the Network Dataset that was just created.\n\n') + print('Instructions:') + print(f'-Expand {gdb}') + print('-Expand the "Colorado" feature dataset') + print('-Right-click the Route_Network dataset that was just created') + print('-Select Properties -> Source Settings -> Group Connectivity') + print('-Change the policy for Route_Connectors from "Honor" to "Override".\n') print('After making this change, uncomment the code below and run it to rebuild the Network to incorporate the change and solve.') - # print('Rebuilding Network...') - # arcpy.na.BuildNetwork(network) - # print('Network rebuilt!') - # - # print('Solving nearest 3 breweries to each fourteener') - # arcpy.na.Solve(closest_facility_lyr) - # print('Routes from each Fourteener to nearest 3 breweries created!') - # if arcpy.CheckExtension('Network') == 'Available': - # arcpy.CheckInExtension('Network') + print('Rebuilding Network...') + arcpy.na.BuildNetwork(network) + print('Network rebuilt!') + + print('Solving nearest 3 breweries to each fourteener') + arcpy.na.Solve(closest_facility_lyr) + print('Routes from each Fourteener to nearest 3 breweries created!') + if arcpy.CheckExtension('Network') == 'Available': + arcpy.CheckInExtension('Network') + + # Summarize the results! + routes_lyr = [lyr for lyr in closest_facility_lyr.listLayers() if lyr.name == 'Routes'][0] + arcpy.management.AddField(in_table=routes_lyr, field_name='Distance_Miles', field_type='DOUBLE') + arcpy.management.CalculateField( + in_table=routes_lyr, + field='Distance_Miles', + expression='!shape.length@meters! * 0.000621371', + expression_type='PYTHON3' + ) + summarize_results(routes_lyr) # endregion +def summarize_results(routes_lyr): + """ Prints out the 3 nearest breweries for each fourteener. """ + + summary = {} + with arcpy.da.SearchCursor(routes_lyr, field_names=['FacilityRank', 'Name', 'Distance_Miles']) as sCursor: + for rank, name, mileage in sCursor: + # Split the 'Name' field into peak and brewery (format "Peak - Brewery") + peak, brewery = name.split(' - ') + + # Initialize the peak entry if it doesn't exist + if peak not in summary: + summary[peak] = {} + + mileage = round(mileage, 1) + + # Store the brewery with its rank + summary[peak][rank] = f'{brewery} ({mileage} miles)' + + for peak, ranked_breweries in summary.items(): + print(peak) + for rank in sorted(ranked_breweries.keys()): + print(f'\t{rank}. {ranked_breweries[rank]}') + + if __name__ == '__main__': main() From cf0d8587711dd09948969f1b10324e9d11e5b5ce Mon Sep 17 00:00:00 2001 From: Ed Conrad Date: Wed, 28 May 2025 18:05:12 -0600 Subject: [PATCH 5/5] added Notebook containing the code --- ...GeoParquet in ArcGIS Pro with DuckDB.ipynb | 896 ++++++++++++++++++ 1 file changed, 896 insertions(+) create mode 100644 geoparquet/How To Read And Query GeoParquet in ArcGIS Pro with DuckDB.ipynb diff --git a/geoparquet/How To Read And Query GeoParquet in ArcGIS Pro with DuckDB.ipynb b/geoparquet/How To Read And Query GeoParquet in ArcGIS Pro with DuckDB.ipynb new file mode 100644 index 0000000..a2bbf9a --- /dev/null +++ b/geoparquet/How To Read And Query GeoParquet in ArcGIS Pro with DuckDB.ipynb @@ -0,0 +1,896 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Python Imports" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import arcpy\n", + "import os\n", + "import duckdb\n", + "from collections import namedtuple" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Environment" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Creating a feature dataset called 'Colorado' in C:\\Users\\EdConrad\\Documents\\ArcGIS\\Projects\\Data_Management\\Data_Management.gdb\n" + ] + } + ], + "source": [ + "pro_project = os.path.join('CURRENT')\n", + "aprx = arcpy.mp.ArcGISProject(pro_project)\n", + "gdb = aprx.defaultGeodatabase\n", + "arcpy.env.workspace = gdb\n", + "feature_datasets = arcpy.ListDatasets()\n", + "wgs84 = arcpy.SpatialReference(4326)\n", + "fd_name = 'Colorado'\n", + "if fd_name not in feature_datasets:\n", + " print(f\"Creating a feature dataset called 'Colorado' in {gdb}\")\n", + " arcpy.management.CreateFeatureDataset(gdb, out_name=fd_name, spatial_reference=wgs84)\n", + "\n", + "out_path = os.path.join(gdb, fd_name)\n", + "arcpy.env.overwriteOutput = True\n", + "\n", + "# Create a DuckDB Database instance\n", + "conn = duckdb.connect()\n", + "conn.sql('install spatial;load spatial;') # Geospatial extension that adds support for working with spatial data and functions\n", + "conn.sql('install httpfs;load httpfs;') # Adds support for reading and writing files over an HTTP(S) or S3 connection\n", + "conn.sql(\"set s3_region='us-west-2';\") # Tells DuckDB which AWS region to use when connecting to S3 buckets\n", + "conn.sql('set enable_object_cache=true;') # Improves performance by caching S3 reads locally for reuse\n", + "\n", + "# Define a BoundingBox named tuple and populate it with Colorado's bounding box values\n", + "BoundingBox = namedtuple(typename='BoundingBox', field_names=['xmin', 'xmax', 'ymin', 'ymax'])\n", + "colorado_bbox = BoundingBox(xmin=-109, xmax=-102, ymin=37, ymax=41)\n", + "# colorado_bbox = BoundingBox(xmin=-108.2, xmax=-104.5, ymin=37, ymax=41)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Get Colorado Fourteeners!" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Creating the Fourteeners feature class...\n", + "\t-Inserted 0 peaks...\n", + "\t-Inserted 5 peaks...\n", + "\t-Inserted 10 peaks...\n", + "\t-Inserted 15 peaks...\n", + "\t-Inserted 20 peaks...\n", + "\t-Inserted 25 peaks...\n", + "\t-Inserted 30 peaks...\n", + "\t-Inserted 35 peaks...\n", + "\t-Inserted 40 peaks...\n", + "\t-Inserted 45 peaks...\n", + "\t-Inserted 50 peaks...\n", + "\t-Inserted 55 peaks...\n", + "\tSummary: 57 points exist in the Fourteeners feature class.\n" + ] + } + ], + "source": [ + "# region Get Colorado Fourteeners!\n", + "# Overture Maps is an open data initiative that provides geospatial data in GeoParquet format,\n", + "# which can be queried efficiently using DuckDB’s built-in Parquet and spatial extensions.\n", + "# Dataset Schema: https://docs.overturemaps.org/schema/reference/base/land/\n", + "# Note 1: GeoParquet stores geometries in a Well-Known Binary (WKB), which means it supports all vector\n", + "# geometry types defined by the OGC Simple Features standard.\n", + "# Note 2: we leverage DuckDB's SQL dialect to select only the columns we want, passing a boundary box to limit the\n", + "# results to Colorado, and we make use of DuckDB's read_parquet() method and the DuckDB spatial extension's\n", + "# ST_AsWKB method to read the WKB\n", + "# Overture Maps regularly updates the data;\n", + "# NOTE: GeoParquet is an immutable file type, meaning new GeoParquet files are created if there's a change\n", + "release = '2025-05-21.0'\n", + "sql = f\"\"\"\n", + " SELECT\n", + " names.primary AS name,\n", + " CAST(elevation AS INT) AS elevation_meters,\n", + " CAST(elevation * 3.28084 AS INT) AS elevation_ft,\n", + " elevation * 3.28084 > 14000 AS Is_Fourteener,\n", + " ST_AsWKB(geometry) AS wkb\n", + " FROM\n", + " read_parquet('s3://overturemaps-us-west-2/release/{release}/theme=base/type=land/*')\n", + " WHERE \n", + " subtype = 'physical' AND class IN ('peak','volcano') AND elevation IS NOT NULL\n", + " AND elevation * 3.28084 > 14000\n", + " AND bbox.xmin BETWEEN {colorado_bbox.xmin} AND {colorado_bbox.xmax}\n", + " AND bbox.ymin BETWEEN {colorado_bbox.ymin} AND {colorado_bbox.ymax}\n", + "\"\"\"\n", + "\n", + "# Get DuckDB Query Relation object\n", + "duck_peaks = conn.sql(sql)\n", + "\n", + "# Insert that data read from GeoParquet into an Esri Feature Class\n", + "import arcpy\n", + "print('Creating the Fourteeners feature class...')\n", + "peaks_fc = arcpy.management.CreateFeatureclass(out_path=out_path,\n", + " out_name='Fourteeners',\n", + " geometry_type='POINT',\n", + " spatial_reference=wgs84).getOutput(0)\n", + "arcpy.management.AddField(in_table=peaks_fc, field_name='Name', field_type='TEXT', field_length=50)\n", + "arcpy.management.AddField(in_table=peaks_fc, field_name='Elevation_meters', field_type='SHORT')\n", + "arcpy.management.AddField(in_table=peaks_fc, field_name='Elevation_ft', field_type='SHORT')\n", + "arcpy.management.AddField(in_table=peaks_fc, field_name='Is_Fourteener', field_type='SHORT')\n", + "with arcpy.da.InsertCursor(peaks_fc, field_names=['Name', 'Elevation_meters', 'Elevation_ft', 'Is_Fourteener',\n", + " 'shape@']) as iCursor:\n", + " # Fetches a single row as a tuple\n", + " raw_row = duck_peaks.fetchone()\n", + " i = 0\n", + " while raw_row:\n", + " if i % 5 == 0:\n", + " print(f'\\t-Inserted {i} peaks...')\n", + "\n", + " # Cast the tuple to list since we need to modify it\n", + " row = list(raw_row)\n", + "\n", + " # Convert the Well-Known Binary (WKB) to an Esri Geometry object\n", + " row[-1] = arcpy.FromWKB(row[-1])\n", + " iCursor.insertRow(row)\n", + " raw_row = duck_peaks.fetchone()\n", + " i += 1\n", + "print(f'\\tSummary: {i} points exist in the Fourteeners feature class.')\n", + "# endregion" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Creating the Breweries feature class...\n", + "\t-Inserted 0 breweries...\n", + "\t-Inserted 50 breweries...\n", + "\t-Inserted 100 breweries...\n", + "\t-Inserted 150 breweries...\n", + "\t-Inserted 200 breweries...\n", + "\t-Inserted 250 breweries...\n", + "\t-Inserted 300 breweries...\n", + "\t-Inserted 350 breweries...\n", + "\t-Inserted 400 breweries...\n", + "\t-Inserted 450 breweries...\n", + "\t-Inserted 500 breweries...\n", + "\t-Inserted 550 breweries...\n", + "\tSummary: 580 points exist in the Breweries feature class.\n" + ] + } + ], + "source": [ + "# region Get Colorado Breweries!\n", + "sql = f\"\"\"\n", + " SELECT\n", + " names.primary AS Name,\n", + " addresses[1].freeform AS Address,\n", + " addresses[1].locality AS City,\n", + " addresses[1].postcode as Zip,\n", + " phones[1] as PhoneNumber,\n", + " ST_AsWKB(geometry) AS wkb\n", + " FROM\n", + " read_parquet('s3://overturemaps-us-west-2/release/{release}/theme=places/type=place/*')\n", + " WHERE \n", + " categories.primary = 'brewery'\n", + " AND bbox.xmin BETWEEN {colorado_bbox.xmin} AND {colorado_bbox.xmax}\n", + " AND bbox.ymin BETWEEN {colorado_bbox.ymin} AND {colorado_bbox.ymax}\n", + " \"\"\"\n", + "\n", + "# Get DuckDB Query Relation object\n", + "duck_breweries = conn.sql(sql)\n", + "\n", + "# Insert that data read from GeoParquet into an Esri Feature Class\n", + "print('Creating the Breweries feature class...')\n", + "breweries_fc = arcpy.management.CreateFeatureclass(out_path=out_path,\n", + " out_name='Breweries',\n", + " geometry_type='POINT',\n", + " spatial_reference=wgs84).getOutput(0)\n", + "arcpy.management.AddField(in_table=breweries_fc, field_name='Name', field_type='TEXT', field_length=100)\n", + "arcpy.management.AddField(in_table=breweries_fc, field_name='Address', field_type='TEXT', field_length=100)\n", + "arcpy.management.AddField(in_table=breweries_fc, field_name='City', field_type='TEXT', field_length=20)\n", + "arcpy.management.AddField(in_table=breweries_fc, field_name='Zip', field_type='TEXT', field_length=10)\n", + "arcpy.management.AddField(in_table=breweries_fc, field_name='Phone_Number', field_type='TEXT', field_length=12)\n", + "\n", + "with arcpy.da.InsertCursor(breweries_fc, field_names=['Name', 'Address', 'City', 'Zip', 'Phone_Number', 'shape@']) as iCursor:\n", + " # Get a single row as a tuple\n", + " raw_row = duck_breweries.fetchone()\n", + " i = 0\n", + " while raw_row:\n", + " if i % 50 == 0:\n", + " print(f'\\t-Inserted {i} breweries...')\n", + "\n", + " # Cast the tuple to list since we need to modify it\n", + " row = list(raw_row)\n", + "\n", + " # Convert the Well-Known Binary (WKB) to an Esri Geometry object\n", + " row[-1] = arcpy.FromWKB(row[-1])\n", + " iCursor.insertRow(row)\n", + " raw_row = duck_breweries.fetchone()\n", + " i += 1\n", + "print(f'\\tSummary: {i} points exist in the Breweries feature class.')\n", + "# endregion" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Creating the Route_Segments feature class...\n", + "\t-Inserted 0 segments...\n", + "\t-Inserted 25,000 segments...\n", + "\t-Inserted 50,000 segments...\n", + "\t-Inserted 75,000 segments...\n", + "\t-Inserted 100,000 segments...\n", + "\t-Inserted 125,000 segments...\n", + "\t-Inserted 150,000 segments...\n", + "\t-Inserted 175,000 segments...\n", + "\t-Inserted 200,000 segments...\n", + "\t-Inserted 225,000 segments...\n", + "\t-Inserted 250,000 segments...\n", + "\t-Inserted 275,000 segments...\n", + "\t-Inserted 300,000 segments...\n", + "\t-Inserted 325,000 segments...\n", + "\t-Inserted 350,000 segments...\n", + "\t-Inserted 375,000 segments...\n", + "\t-Inserted 400,000 segments...\n", + "\t-Inserted 425,000 segments...\n", + "\t-Inserted 450,000 segments...\n", + "\t-Inserted 475,000 segments...\n", + "\t-Inserted 500,000 segments...\n", + "\t-Inserted 525,000 segments...\n", + "\t-Inserted 550,000 segments...\n", + "\t-Inserted 575,000 segments...\n", + "\t-Inserted 600,000 segments...\n", + "\t-Inserted 625,000 segments...\n", + "\t-Inserted 650,000 segments...\n", + "\t-Inserted 675,000 segments...\n", + "\t-Inserted 700,000 segments...\n", + "\t-Inserted 725,000 segments...\n", + "\t-Inserted 750,000 segments...\n", + "\t-Inserted 775,000 segments...\n", + "\t-Inserted 800,000 segments...\n", + "\t-Inserted 825,000 segments...\n", + "\t-Inserted 850,000 segments...\n", + "\t-Inserted 875,000 segments...\n", + "\t-Inserted 900,000 segments...\n", + "\t-Inserted 925,000 segments...\n", + "\t-Inserted 950,000 segments...\n", + "\t-Inserted 975,000 segments...\n", + "\t-Inserted 1,000,000 segments...\n", + "\t-Inserted 1,025,000 segments...\n", + "\t-Inserted 1,050,000 segments...\n", + "\t-Inserted 1,075,000 segments...\n", + "\t-Inserted 1,100,000 segments...\n", + "\t-Inserted 1,125,000 segments...\n", + "\t-Inserted 1,150,000 segments...\n", + "\t-Inserted 1,175,000 segments...\n", + "\tSummary: 1,175,907 polylines exist in the Route_Segments feature class.\n", + "Creating the Route_Connectors feature class...\n", + "\t-Inserted 0 connectors...\n", + "\t-Inserted 25,000 connectors...\n", + "\t-Inserted 50,000 connectors...\n", + "\t-Inserted 75,000 connectors...\n", + "\t-Inserted 100,000 connectors...\n", + "\t-Inserted 125,000 connectors...\n", + "\t-Inserted 150,000 connectors...\n", + "\t-Inserted 175,000 connectors...\n", + "\t-Inserted 200,000 connectors...\n", + "\t-Inserted 225,000 connectors...\n", + "\t-Inserted 250,000 connectors...\n", + "\t-Inserted 275,000 connectors...\n", + "\t-Inserted 300,000 connectors...\n", + "\t-Inserted 325,000 connectors...\n", + "\t-Inserted 350,000 connectors...\n", + "\t-Inserted 375,000 connectors...\n", + "\t-Inserted 400,000 connectors...\n", + "\t-Inserted 425,000 connectors...\n", + "\t-Inserted 450,000 connectors...\n", + "\t-Inserted 475,000 connectors...\n", + "\t-Inserted 500,000 connectors...\n", + "\t-Inserted 525,000 connectors...\n", + "\t-Inserted 550,000 connectors...\n", + "\t-Inserted 575,000 connectors...\n", + "\t-Inserted 600,000 connectors...\n", + "\t-Inserted 625,000 connectors...\n", + "\t-Inserted 650,000 connectors...\n", + "\t-Inserted 675,000 connectors...\n", + "\t-Inserted 700,000 connectors...\n", + "\t-Inserted 725,000 connectors...\n", + "\t-Inserted 750,000 connectors...\n", + "\t-Inserted 775,000 connectors...\n", + "\t-Inserted 800,000 connectors...\n", + "\t-Inserted 825,000 connectors...\n", + "\t-Inserted 850,000 connectors...\n", + "\t-Inserted 875,000 connectors...\n", + "\t-Inserted 900,000 connectors...\n", + "\t-Inserted 925,000 connectors...\n", + "\t-Inserted 950,000 connectors...\n", + "\t-Inserted 975,000 connectors...\n", + "\t-Inserted 1,000,000 connectors...\n", + "\t-Inserted 1,025,000 connectors...\n", + "\t-Inserted 1,050,000 connectors...\n", + "\t-Inserted 1,075,000 connectors...\n", + "\t-Inserted 1,100,000 connectors...\n", + "\t-Inserted 1,125,000 connectors...\n", + "\t-Inserted 1,150,000 connectors...\n", + "\t-Inserted 1,175,000 connectors...\n", + "\t-Inserted 1,200,000 connectors...\n", + "\t-Inserted 1,225,000 connectors...\n", + "\t-Inserted 1,250,000 connectors...\n", + "\t-Inserted 1,275,000 connectors...\n", + "\t-Inserted 1,300,000 connectors...\n", + "\t-Inserted 1,325,000 connectors...\n", + "\t-Inserted 1,350,000 connectors...\n", + "\t-Inserted 1,375,000 connectors...\n", + "\t-Inserted 1,400,000 connectors...\n", + "\t-Inserted 1,425,000 connectors...\n", + "\t-Inserted 1,450,000 connectors...\n", + "\t-Inserted 1,475,000 connectors...\n", + "\t-Inserted 1,500,000 connectors...\n", + "\t-Inserted 1,525,000 connectors...\n", + "\t-Inserted 1,550,000 connectors...\n", + "\t-Inserted 1,575,000 connectors...\n", + "\t-Inserted 1,600,000 connectors...\n", + "\t-Inserted 1,625,000 connectors...\n", + "\t-Inserted 1,650,000 connectors...\n", + "\t-Inserted 1,675,000 connectors...\n", + "\t-Inserted 1,700,000 connectors...\n", + "\t-Inserted 1,725,000 connectors...\n", + "\tSummary: 1,730,787 points exist in the Route_Connectors feature class.\n" + ] + } + ], + "source": [ + "# region Get Colorado Road Segments!\n", + "# According to the docs, the definition for subtype of 'road':\n", + "# \"A road segment represents a section of any kind of road, street or path, including a dedicated path for walking or cycling, but excluding a railway.\"\n", + "# As of 5/28/25 the names.primary is returning segment type - which differs from the docs.\n", + "sql = f\"\"\"\n", + " SELECT\n", + " names.primary as SegmentType, \n", + " ST_AsWKB(geometry) AS wkb\n", + " FROM\n", + " read_parquet('s3://overturemaps-us-west-2/release/{release}/theme=transportation/type=segment/*')\n", + " WHERE \n", + " subtype = 'road'\n", + " AND bbox.xmin BETWEEN {colorado_bbox.xmin} AND {colorado_bbox.xmax}\n", + " AND bbox.ymin BETWEEN {colorado_bbox.ymin} AND {colorado_bbox.ymax}\n", + " \"\"\"\n", + "\n", + "# Get DuckDB Query Relation object\n", + "duck_segments = conn.sql(sql)\n", + "\n", + "# Insert that data read from GeoParquet into an Esri Feature Class\n", + "print('Creating the Route_Segments feature class...')\n", + "segments_fc = arcpy.management.CreateFeatureclass(out_path=out_path,\n", + " out_name='Route_Segments',\n", + " geometry_type='POLYLINE',\n", + " spatial_reference=wgs84).getOutput(0)\n", + "arcpy.management.AddField(in_table=segments_fc, field_name='Segment_Type', field_type='TEXT', field_length=20)\n", + "with arcpy.da.InsertCursor(segments_fc, field_names=['shape@', 'Segment_Type']) as iCursor:\n", + " # Get a single row as a tuple\n", + " raw_row = duck_segments.fetchone()\n", + " i = 0\n", + " while raw_row:\n", + " if i % 25_000 == 0:\n", + " print(f'\\t-Inserted {i:,} segments...')\n", + "\n", + " # Cast the tuple to list since we need to modify it\n", + " row = list(raw_row)\n", + "\n", + " # Convert the Well-Known Binary (WKB) string to an Esri Geometry object\n", + " row[0] = arcpy.FromWKB(row[-1])\n", + " iCursor.insertRow(row)\n", + " raw_row = duck_segments.fetchone()\n", + " i += 1\n", + "print(f'\\tSummary: {i:,} polylines exist in the Route_Segments feature class.')\n", + "\n", + "# region Get Colorado Road Connectors!\n", + "sql = f\"\"\"\n", + " SELECT\n", + " ST_AsWKB(geometry) AS wkb\n", + " FROM\n", + " read_parquet('s3://overturemaps-us-west-2/release/{release}/theme=transportation/type=connector/*')\n", + " WHERE \n", + " bbox.xmin BETWEEN {colorado_bbox.xmin} AND {colorado_bbox.xmax}\n", + " AND bbox.ymin BETWEEN {colorado_bbox.ymin} AND {colorado_bbox.ymax}\n", + " \"\"\"\n", + "\n", + "# Get DuckDB Query Relation object\n", + "duck_connectors = conn.sql(sql)\n", + "\n", + "# Insert that data read from GeoParquet into an Esri Feature Class\n", + "print('Creating the Route_Connectors feature class...')\n", + "connectors_fc = arcpy.management.CreateFeatureclass(out_path=out_path,\n", + " out_name='Route_Connectors',\n", + " geometry_type='POINT',\n", + " spatial_reference=wgs84).getOutput(0)\n", + "\n", + "with arcpy.da.InsertCursor(connectors_fc, field_names=['shape@']) as iCursor:\n", + " # Get a single row as a tuple\n", + " raw_row = duck_connectors.fetchone()\n", + " i = 0\n", + " while raw_row:\n", + " if i % 25_000 == 0:\n", + " print(f'\\t-Inserted {i:,} connectors...')\n", + "\n", + " # Cast the tuple to list since we need to modify it\n", + " row = list(raw_row)\n", + "\n", + " # Convert the Well-Known Binary (WKB) string to an Esri Geometry object\n", + " row[0] = arcpy.FromWKB(row[-1])\n", + " iCursor.insertRow(row)\n", + " raw_row = duck_connectors.fetchone()\n", + " i += 1\n", + "print(f'\\tSummary: {i:,} points exist in the Route_Connectors feature class.')\n", + "# endregion" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Close DuckDB" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "conn.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Network Analyst - Closest Facility Analysis" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Creating Route_Network...\n", + "Creating Closest Facility Analysis Layers\n", + "Building Network...\n", + "Network built!\n", + "Adding Breweries as the Facilities\n", + "Adding Fourteeners as the Incidents\n", + "\n", + "\n", + "Now you need to make a manual property change to the Network Dataset that was just created.\n", + "\n", + "\n", + "Instructions:\n", + "-Expand C:\\Users\\EdConrad\\Documents\\ArcGIS\\Projects\\Data_Management\\Data_Management.gdb\n", + "-Expand the \"Colorado\" feature dataset\n", + "-Right-click the Route_Network dataset that was just created\n", + "-Select Properties -> Source Settings -> Group Connectivity\n", + "-Change the policy for Route_Connectors from \"Honor\" to \"Override\".\n", + "\n", + "After making this change, uncomment the code below and run it to rebuild the Network to incorporate the change and solve.\n" + ] + } + ], + "source": [ + "# region Network Analyst section\n", + "# Check out Network Analyst license if available. Fail if the Network Analyst license is not available.\n", + "if arcpy.CheckExtension('Network') == 'Available':\n", + " arcpy.CheckOutExtension('Network')\n", + " import arcpy.na\n", + "else:\n", + " raise arcpy.ExecuteError('Network Analyst Extension license is not available.')\n", + "\n", + "print('Creating Route_Network...')\n", + "\n", + "out_path = \"C:\\\\Users\\\\EdConrad\\\\Documents\\\\ArcGIS\\\\Projects\\\\Data_Management\\\\Testing.gdb\\\\Colorado\"\n", + "\n", + "arcpy.na.CreateNetworkDataset(feature_dataset=out_path,\n", + " out_name='Route_Network',\n", + " source_feature_class_names='Route_Segments;Route_Connectors',\n", + " elevation_model='NO_ELEVATION')\n", + "network = os.path.join(out_path, 'Route_Network')\n", + "\n", + "print('Creating Closest Facility Analysis Layers')\n", + "closest_facility_lyr = arcpy.na.MakeClosestFacilityAnalysisLayer(network_data_source=network,\n", + " layer_name='ClosestBreweriesToFourteeners',\n", + " travel_direction='TO_FACILITIES',\n", + " number_of_facilities_to_find=3).getOutput(0)\n", + "print('Building Network...')\n", + "arcpy.na.BuildNetwork(network)\n", + "print('Network built!')\n", + "\n", + "# Add Breweries as the \"Facilities\"\n", + "print('Adding Breweries as the Facilities')\n", + "arcpy.na.AddLocations(\n", + " in_network_analysis_layer=closest_facility_lyr,\n", + " sub_layer='Facilities',\n", + " in_table=breweries_fc,\n", + " search_tolerance='100 Meters', # NOTE: value determined after some trial-and-error\n", + " match_type='MATCH_TO_CLOSEST',\n", + " append='APPEND',\n", + " snap_to_position_along_network='SNAP',\n", + " snap_offset='0 Meters'\n", + ")\n", + "\n", + "# Add Fourteeners as the \"Incidents\"\n", + "print('Adding Fourteeners as the Incidents')\n", + "arcpy.na.AddLocations(\n", + " in_network_analysis_layer=closest_facility_lyr,\n", + " sub_layer='Incidents',\n", + " in_table=peaks_fc,\n", + " search_tolerance='1500 Meters', # NOTE: value determined after some trial-and-error\n", + " match_type='MATCH_TO_CLOSEST',\n", + " append='APPEND',\n", + " snap_to_position_along_network='SNAP',\n", + " snap_offset='0 Meters'\n", + ")\n", + "\n", + "print('\\n\\nNow you need to make a manual property change to the Network Dataset that was just created.\\n\\n')\n", + "print('Instructions:')\n", + "print(f'-Expand {gdb}')\n", + "print('-Expand the \"Colorado\" feature dataset')\n", + "print('-Right-click the Route_Network dataset that was just created')\n", + "print('-Select Properties -> Source Settings -> Group Connectivity')\n", + "print('-Change the policy for Route_Connectors from \"Honor\" to \"Override\".\\n')\n", + "print('After making this change, uncomment the code below and run it to rebuild the Network to incorporate the change and solve.')\n", + "# endregion" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## After making the manual change noted above, uncomment lines below and run" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "El Diente Peak\n", + "\t1. Telluride Brewing Co. (24.7 miles)\n", + "\t2. Avalanche Brewing Co (34.7 miles)\n", + "\t3. Golden Block Brewery (34.8 miles)\n", + "Mount Wilson\n", + "\t1. Telluride Brewing Co. (24.9 miles)\n", + "\t2. Avalanche Brewing Co (35.0 miles)\n", + "\t3. Golden Block Brewery (35.0 miles)\n", + "Wilson Peak\n", + "\t1. Telluride Brewing Co. (15.1 miles)\n", + "\t2. Avalanche Brewing Co (28.3 miles)\n", + "\t3. Golden Block Brewery (28.3 miles)\n", + "Mount Eolus\n", + "\t1. Golden Block Brewery (39.3 miles)\n", + "\t2. Avalanche Brewing Co (39.4 miles)\n", + "\t3. Bottom Shelf Brewery (42.2 miles)\n", + "North Eolus\n", + "\t1. Golden Block Brewery (39.2 miles)\n", + "\t2. Avalanche Brewing Co (39.3 miles)\n", + "\t3. Bottom Shelf Brewery (42.0 miles)\n", + "Windom Peak\n", + "\t1. Golden Block Brewery (39.2 miles)\n", + "\t2. Avalanche Brewing Co (39.3 miles)\n", + "\t3. Bottom Shelf Brewery (42.0 miles)\n", + "Sunlight Peak\n", + "\t1. Golden Block Brewery (39.2 miles)\n", + "\t2. Avalanche Brewing Co (39.3 miles)\n", + "\t3. Bottom Shelf Brewery (42.0 miles)\n", + "Handies Peak\n", + "\t1. Golden Block Brewery (15.7 miles)\n", + "\t2. Avalanche Brewing Co (15.8 miles)\n", + "\t3. Ouray Brewery (18.0 miles)\n", + "Sunshine Peak\n", + "\t1. Lake City Brewing Company (22.6 miles)\n", + "\t2. Golden Block Brewery (23.6 miles)\n", + "\t3. Avalanche Brewing Co (23.6 miles)\n", + "Redcloud Peak\n", + "\t1. Lake City Brewing Company (23.2 miles)\n", + "\t2. Golden Block Brewery (24.2 miles)\n", + "\t3. Avalanche Brewing Co (24.2 miles)\n", + "Mount Sneffels\n", + "\t1. Ouray Brewery (10.7 miles)\n", + "\t2. Mr. Grumpy Pants Brewing (10.8 miles)\n", + "\t3. Telluride Brewing Co. (14.7 miles)\n", + "Wetterhorn Peak\n", + "\t1. Lake City Brewing Company (14.7 miles)\n", + "\t2. Ouray Brewery (18.6 miles)\n", + "\t3. Mr. Grumpy Pants Brewing (18.7 miles)\n", + "Uncompahgre Peak\n", + "\t1. Lake City Brewing Company (12.8 miles)\n", + "\t2. Ouray Brewery (21.7 miles)\n", + "\t3. Mr. Grumpy Pants Brewing (21.8 miles)\n", + "San Luis Peak\n", + "\t1. Spare Keg Brewerks- Creede (14.3 miles)\n", + "\t2. Lake City Brewing Company (35.0 miles)\n", + "\t3. Three Barrel Brewing Co. (51.8 miles)\n", + "Snowmass Mountain\n", + "\t1. The Eldo Brewpub & Venue (23.5 miles)\n", + "\t2. Irwin Brewing Company (23.5 miles)\n", + "\t3. Aspen Brewing Company (24.0 miles)\n", + "Capitol Peak\n", + "\t1. Aspen Brewing Company (21.6 miles)\n", + "\t2. Capitol Creek Brewery (23.1 miles)\n", + "\t3. Aspen Brewing Company (23.8 miles)\n", + "Maroon Peak\n", + "\t1. Aspen Brewing Company (16.0 miles)\n", + "\t2. Aspen Brewing Company (16.4 miles)\n", + "\t3. The Eldo Brewpub & Venue (23.6 miles)\n", + "Pyramid Peak\n", + "\t1. Aspen Brewing Company (13.3 miles)\n", + "\t2. Aspen Brewing Company (13.8 miles)\n", + "\t3. The Eldo Brewpub & Venue (25.0 miles)\n", + "Little Bear Peak\n", + "\t1. Spare Keg Brewerks (26.6 miles)\n", + "\t2. San Luis Valley Brewing Co. (26.6 miles)\n", + "\t3. The Colorado Farm Brewery (37.2 miles)\n", + "Ellingwood Point\n", + "\t1. Spare Keg Brewerks (28.0 miles)\n", + "\t2. San Luis Valley Brewing Co. (28.0 miles)\n", + "\t3. The Colorado Farm Brewery (38.6 miles)\n", + "Blanca Peak\n", + "\t1. Spare Keg Brewerks (27.8 miles)\n", + "\t2. San Luis Valley Brewing Co. (27.8 miles)\n", + "\t3. The Colorado Farm Brewery (38.4 miles)\n", + "Crestone Peak\n", + "\t1. Spare Keg Brewerks (47.2 miles)\n", + "\t2. San Luis Valley Brewing Co. (47.2 miles)\n", + "\t3. C&S Brewery & Bistro (50.0 miles)\n", + "Crestone Needle\n", + "\t1. Spare Keg Brewerks (47.6 miles)\n", + "\t2. San Luis Valley Brewing Co. (47.6 miles)\n", + "\t3. C&S Brewery & Bistro (49.0 miles)\n", + "Culebra Peak\n", + "\t1. Spare Keg Brewerks (57.6 miles)\n", + "\t2. San Luis Valley Brewing Co. (57.6 miles)\n", + "\t3. The Colorado Farm Brewery (61.6 miles)\n", + "Mount Lindsey\n", + "\t1. Mountain Merman Brewing Company (45.7 miles)\n", + "\t2. Crafty Canary Brewery (53.7 miles)\n", + "\t3. Spare Keg Brewerks (57.3 miles)\n", + "Tabeguache Peak\n", + "\t1. Elevation Beer Company (15.0 miles)\n", + "\t2. Soulcraft Beer (17.7 miles)\n", + "\t3. Tres Litros Beer Company (18.1 miles)\n", + "Mount Shavano\n", + "\t1. Elevation Beer Company (14.1 miles)\n", + "\t2. Soulcraft Beer (16.8 miles)\n", + "\t3. Tres Litros Beer Company (17.2 miles)\n", + "Kit Carson Mountain\n", + "\t1. Spare Keg Brewerks (51.9 miles)\n", + "\t2. San Luis Valley Brewing Co. (52.0 miles)\n", + "\t3. Coors Adolph Co (59.3 miles)\n", + "Challenger Point\n", + "\t1. Spare Keg Brewerks (51.4 miles)\n", + "\t2. San Luis Valley Brewing Co. (51.4 miles)\n", + "\t3. Coors Adolph Co (58.7 miles)\n", + "Humboldt Peak\n", + "\t1. C&S Brewery & Bistro (49.8 miles)\n", + "\t2. Spare Keg Brewerks (49.9 miles)\n", + "\t3. San Luis Valley Brewing Co. (49.9 miles)\n", + "Castle Peak\n", + "\t1. Aspen Brewing Company (19.7 miles)\n", + "\t2. Aspen Brewing Company (21.2 miles)\n", + "\t3. Irwin Brewing Company (23.4 miles)\n", + "Conundrum Peak\n", + "\t1. Aspen Brewing Company (19.8 miles)\n", + "\t2. Aspen Brewing Company (21.2 miles)\n", + "\t3. Irwin Brewing Company (23.4 miles)\n", + "Mount Antero\n", + "\t1. Eddyline Brewery (23.0 miles)\n", + "\t2. Elevation Beer Company (23.8 miles)\n", + "\t3. Eddyline Restaurant at Southmain (23.8 miles)\n", + "Mount Princeton\n", + "\t1. Eddyline Brewery (14.4 miles)\n", + "\t2. Eddyline Restaurant at Southmain (15.2 miles)\n", + "\t3. Elevation Beer Company (26.8 miles)\n", + "Mount Yale\n", + "\t1. Eddyline Brewery (12.2 miles)\n", + "\t2. Eddyline Restaurant at Southmain (13.0 miles)\n", + "\t3. Elevation Beer Company (36.4 miles)\n", + "Huron Peak\n", + "\t1. Cloud City Modern Mead (31.6 miles)\n", + "\t2. Eddyline Brewery (32.1 miles)\n", + "\t3. Two Mile Brewing Co. (32.2 miles)\n", + "Missouri Mountain\n", + "\t1. Eddyline Brewery (23.5 miles)\n", + "\t2. Eddyline Restaurant at Southmain (24.3 miles)\n", + "\t3. Cloud City Modern Mead (29.8 miles)\n", + "Mount Belford\n", + "\t1. Eddyline Brewery (23.0 miles)\n", + "\t2. Eddyline Restaurant at Southmain (23.8 miles)\n", + "\t3. Cloud City Modern Mead (29.5 miles)\n", + "Mount Harvard\n", + "\t1. Eddyline Brewery (14.9 miles)\n", + "\t2. Eddyline Restaurant at Southmain (15.7 miles)\n", + "\t3. Elevation Beer Company (39.2 miles)\n", + "Mount Columbia\n", + "\t1. Eddyline Brewery (13.1 miles)\n", + "\t2. Eddyline Restaurant at Southmain (13.9 miles)\n", + "\t3. Elevation Beer Company (37.4 miles)\n", + "Mount Oxford\n", + "\t1. Eddyline Brewery (21.6 miles)\n", + "\t2. Eddyline Restaurant at Southmain (22.4 miles)\n", + "\t3. Cloud City Modern Mead (28.2 miles)\n", + "La Plata Peak\n", + "\t1. Cloud City Modern Mead (28.7 miles)\n", + "\t2. Two Mile Brewing Co. (29.3 miles)\n", + "\t3. Aspen Brewing Company (31.0 miles)\n", + "Mount Elbert\n", + "\t1. Cloud City Modern Mead (15.1 miles)\n", + "\t2. Two Mile Brewing Co. (15.6 miles)\n", + "\t3. Eddyline Brewery (31.3 miles)\n", + "Mount Massive\n", + "\t1. Cloud City Modern Mead (14.2 miles)\n", + "\t2. Two Mile Brewing Co. (14.7 miles)\n", + "\t3. South Park Lager Beer Brewery (35.6 miles)\n", + "Mount Sherman\n", + "\t1. Two Mile Brewing Co. (9.5 miles)\n", + "\t2. Cloud City Modern Mead (9.7 miles)\n", + "\t3. South Park Lager Beer Brewery (14.4 miles)\n", + "Mount Democrat\n", + "\t1. South Park Lager Beer Brewery (12.6 miles)\n", + "\t2. Highside Brewing (13.4 miles)\n", + "\t3. South Park Brewing Colorado (13.5 miles)\n", + "Mount Bross\n", + "\t1. South Park Lager Beer Brewery (12.4 miles)\n", + "\t2. Highside Brewing (13.2 miles)\n", + "\t3. South Park Brewing Colorado (13.2 miles)\n", + "Mount Cameron\n", + "\t1. South Park Lager Beer Brewery (12.9 miles)\n", + "\t2. Highside Brewing (13.7 miles)\n", + "\t3. South Park Brewing Colorado (13.8 miles)\n", + "Mount Lincoln\n", + "\t1. South Park Lager Beer Brewery (13.5 miles)\n", + "\t2. Highside Brewing (14.3 miles)\n", + "\t3. South Park Brewing Colorado (14.3 miles)\n", + "Pikes Peak\n", + "\t1. Manitou Brewing Company (11.5 miles)\n", + "\t2. WestFax Springs (14.0 miles)\n", + "\t3. Paradox Beer Company (14.2 miles)\n", + "Mount of the Holy Cross\n", + "\t1. Vail Brewing Company (19.3 miles)\n", + "\t2. Vail Brewing Company (22.8 miles)\n", + "\t3. Highside Brewing (40.6 miles)\n", + "Quandary Peak\n", + "\t1. Breckenridge Brewery (10.9 miles)\n", + "\t2. Broken Compass Brewery (11.0 miles)\n", + "\t3. Breckenridge Brew Pub (11.2 miles)\n", + "Torreys Peak\n", + "\t1. Guanella Pass Brewing Company (12.6 miles)\n", + "\t2. Steep Brewing & Coffee Company (12.9 miles)\n", + "\t3. Cabin Creek Brewing (13.3 miles)\n", + "Grays Peak\n", + "\t1. Steep Brewing & Coffee Company (12.3 miles)\n", + "\t2. Guanella Pass Brewing Company (12.8 miles)\n", + "\t3. Cabin Creek Brewing (13.9 miles)\n", + "Mount Bierstadt\n", + "\t1. Guanella Pass Brewing Company (14.2 miles)\n", + "\t2. Cabin Creek Brewing (15.7 miles)\n", + "\t3. Tommyknocker Brewery & Pub (18.3 miles)\n", + "Mount Blue Sky\n", + "\t1. Guanella Pass Brewing Company (15.1 miles)\n", + "\t2. Cabin Creek Brewing (16.6 miles)\n", + "\t3. Tommyknocker Brewery & Pub (17.0 miles)\n", + "Longs Peak\n", + "\t1. Rock Cut Brewing Company (15.0 miles)\n", + "\t2. Estes Park Brewery (15.0 miles)\n", + "\t3. Lumpy Ridge Brewing Company (15.7 miles)\n" + ] + } + ], + "source": [ + "# print('Rebuilding Network...')\n", + "# arcpy.na.BuildNetwork(network)\n", + "# print('Network rebuilt!')\n", + "\n", + "# print('Solving nearest 3 breweries to each fourteener')\n", + "# arcpy.na.Solve(closest_facility_lyr)\n", + "# print('Routes from each Fourteener to nearest 3 breweries created!')\n", + "# if arcpy.CheckExtension('Network') == 'Available':\n", + "# arcpy.CheckInExtension('Network')\n", + "\n", + "# def summarize_results(routes_lyr):\n", + "# \"\"\" Prints out the 3 nearest breweries for each fourteener. \"\"\"\n", + "\n", + "# summary = {}\n", + "# with arcpy.da.SearchCursor(routes_lyr, field_names=['FacilityRank', 'Name', 'Distance_Miles']) as sCursor:\n", + "# for rank, name, mileage in sCursor:\n", + "# # Split the 'Name' field into peak and brewery (format \"Peak - Brewery\")\n", + "# peak, brewery = name.split(' - ')\n", + "\n", + "# # Initialize the peak entry if it doesn't exist\n", + "# if peak not in summary:\n", + "# summary[peak] = {}\n", + "\n", + "# mileage = round(mileage, 1)\n", + " \n", + "# # Store the brewery with its rank\n", + "# summary[peak][rank] = f'{brewery} ({mileage} miles)'\n", + "\n", + "# for peak, ranked_breweries in summary.items():\n", + "# print(peak)\n", + "# for rank in sorted(ranked_breweries.keys()):\n", + "# print(f'\\t{rank}. {ranked_breweries[rank]}')\n", + "\n", + "\n", + "# routes_lyr = [lyr for lyr in closest_facility_lyr.listLayers() if lyr.name == 'Routes'][0]\n", + "# arcpy.management.AddField(in_table=routes_lyr, field_name='Distance_Miles', field_type='DOUBLE')\n", + "# arcpy.management.CalculateField(\n", + "# in_table=routes_lyr,\n", + "# field='Distance_Miles',\n", + "# expression='!shape.length@meters! * 0.000621371',\n", + "# expression_type='PYTHON3'\n", + "# )\n", + "\n", + "# # Summarize the results!\n", + "# summarize_results(routes_lyr)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ArcGISPro", + "language": "python", + "name": "python3" + }, + "language_info": { + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "version": "3.11.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}