SQL Server holds spatial data in geography and geometry columns. This post exports two layers of the BD TOPO, IGN's topographic database of France, from SQL Server to Parquet with FastBCP's adbc_mssql connection, and queries the files with DuckDB and GeoPandas. The spatial columns come out as GeoParquet: the geometries are written as WKB, and a geo block in the file footer records which column is spatial, how it is encoded, and in which Coordinate Reference System (CRS). No conversion is needed in the query, and none after the export. FastBCP currently produces GeoParquet 1.1, not 2.0.
The dataset
We use the BD TOPO delivery for the Pyrénées-Orientales department (D066), a 1.3 GB GeoPackage in Lambert-93 (EPSG:2154) with 59 layers, and keep two of them:
troncon_de_route: 210,250 road segments (3D LineStrings)batiment: 451,646 buildings (3D MultiPolygons)
The tables in SQL Server
Both layers already live in SQL Server 2022, running in Docker on the same machine: roads in a geography column, buildings in a geometry column. They were loaded beforehand by a short Python script, which is not the subject of this post. That script flattened the geometries to 2D on the way in: BD TOPO ships 3D coordinates, but Z was discarded in the loading script for simplicity reasons. The export path itself preserves Z and M when the source carries them. M is a fourth ordinate, a measure attached to each vertex, usually a distance along a route so that a position can be given as a kilometer point rather than as coordinates.
The two column types are not interchangeable, and the choice follows from the CRS. geography only accepts a geodetic SRID (Spatial Reference Identifier: the integer code identifying the CRS). The catalog view sys.spatial_reference_systems, present in every database, lists the 393 SRIDs a geography column will accept:
SELECT COUNT(*) AS supported,
SUM(IIF(spatial_reference_id = 4326, 1, 0)) AS has_4326,
SUM(IIF(spatial_reference_id = 2154, 1, 0)) AS has_2154
FROM sys.spatial_reference_systems;supported has_4326 has_2154
--------- -------- --------
393 1 0Lambert-93 (EPSG:2154) is not among them, and cannot be. The dividing line is the projection: a geographic CRS gives a position as longitude and latitude, angles on an ellipsoid, while a projected CRS such as Lambert-93 flattens that ellipsoid onto a plane and gives the position in meters. Lambert-93 is derived from a geographic CRS, but the projection step is what puts it out of reach: geography measures on the ellipsoid, in great-circle distances and geodesic areas, and needs angular coordinates. Every SRID in the sys.spatial_reference_systems is an unprojected system.
The roads were reprojected from Lambert-93 to WGS 84 (EPSG:4326) on the way in, to fit in a geography column. The buildings keep the source's native Lambert-93, in meters, and go into a geometry column, which places no constraint on the SRID and treats coordinates as planar. The export below covers both column types and both kinds of CRS.
The queries below are plain T-SQL and run from any client. Each table has an int IDENTITY key, a few nvarchar and float attributes, and one spatial column:
SELECT COLUMN_NAME, DATA_TYPE
FROM INFORMATION_SCHEMA.COLUMNS
WHERE TABLE_NAME = 'troncon_de_route'
ORDER BY ORDINAL_POSITION;COLUMN_NAME DATA_TYPE
-------------------- ---------
id int
cleabs nvarchar
nature nvarchar
importance nvarchar
nombre_de_voies float
largeur_de_chaussee float
vitesse_moyenne_vl int
sens_de_circulation nvarchar
insee_commune_gauche nvarchar
insee_commune_droite nvarchar
geom geographybatiment has the same overall shape, with its own attributes and a geom column typed geometry instead of geography:
SELECT COLUMN_NAME, DATA_TYPE
FROM INFORMATION_SCHEMA.COLUMNS
WHERE TABLE_NAME = 'batiment'
ORDER BY ORDINAL_POSITION;COLUMN_NAME DATA_TYPE
-------------------- ---------
id int
cleabs nvarchar
nature nvarchar
usage_1 nvarchar
hauteur float
nombre_d_etages float
nombre_de_logements float
materiaux_des_murs nvarchar
geom geometryThe id identity column gives FastBCP an obvious distribution key for the parallel partitioning used below.
Row counts and the SRID carried by each geo column:
SELECT 'troncon_de_route' AS tbl, COUNT(*) AS rows, MIN(geom.STSrid) AS srid
FROM dbo.troncon_de_route
UNION ALL
SELECT 'batiment', COUNT(*), MIN(geom.STSrid) FROM dbo.batiment;tbl rows srid
---------------- ------ ----
troncon_de_route 210250 4326
batiment 451646 2154Parallel export with FastBCP
FastBCP exports each table with 8 parallel readers, partitioned on the id column with the RangeId method. RangeId reads only the key's minimum and maximum, then hands each reader an equal-width id range that rides the identity column. Because id is a dense IDENTITY(1,1), the eight ranges hold near-equal row counts:
FastBCP -C adbc_mssql -S localhost -U francois -X '***' \
-I bdtopo -s dbo -T troncon_de_route \
-m RangeId -c id -p 8 \
-D parquet/roads --fileoutput "{sourcetable}.parquet"
FastBCP -C adbc_mssql -S localhost -U francois -X '***' \
-I bdtopo -s dbo -T batiment \
-m RangeId -c id -p 8 \
-D parquet/buildings --fileoutput "{sourcetable}.parquet"The flags, in the order they appear (the full list is in the CLI reference):
-C adbc_mssql: connection type, here the ADBC SQL Server driver.-S,-U,-X: server, user, password.-I bdtopo: source database.-s dboand-T troncon_de_route: source schema and table. A query can be passed with-qinstead.-m RangeId: parallelism method. Others includeNtile,DataDriven,PhyslocandRandom; the default isNone, meaning a single reader.-c id: the column used to distribute the data across readers.RangeIdrequires it to be numeric.-p 8: degree of parallelism, so eight concurrent readers.-D parquet/roads: output directory, local or on object storage.--fileoutput "{sourcetable}.parquet": output file name, where{sourcetable}expands to the table name. The extension selects the format: csv, tsv, json, bson, parquet, xlsx or binary.
Each table lands as eight Parquet files, one per RangeId partition, Zstd-compressed level 3 by default:
$ ls parquet/roads/
troncon_de_route_chunk_000.parquet
troncon_de_route_chunk_001.parquet
troncon_de_route_chunk_002.parquet
troncon_de_route_chunk_003.parquet
troncon_de_route_chunk_004.parquet
troncon_de_route_chunk_005.parquet
troncon_de_route_chunk_006.parquet
troncon_de_route_chunk_007.parquetBoth runs finish in well under a second:
troncon_de_route: 210250 rows - Total time : Elapsed=538 ms - 390,717 rows/s
batiment: 451646 rows - Total time : Elapsed=635 ms - 710,774 rows/sThe two tables, 661,896 rows and as many geometries, are on disk in a little over 1.1 s. The 16 output files weigh 58 MB, against 1.3 GB for the source GeoPackage, which also holds 57 other layers:
parquet/roads: 8 files, 38 MB total
parquet/buildings: 8 files, 20 MB totalQuerying the Parquet files with DuckDB
DuckDB globs the partition files and reads the GeoParquet metadata with its spatial extension. The geo block in the footer is enough for geom to arrive as a GEOMETRY with its CRS attached, GEOMETRY('OGC:CRS84') for the roads and GEOMETRY('EPSG:2154') for the buildings, so no ST_GeomFromWKB call appears below.
Road lengths in meters require reprojecting the WGS 84 coordinates back to Lambert-93 before ST_Length:
INSTALL spatial; LOAD spatial;
SELECT nature,
count(*) AS segments,
round(sum(ST_Length(ST_Transform(geom,
'EPSG:4326', 'EPSG:2154', always_xy := true))) / 1000) AS km
FROM read_parquet('parquet/roads/*.parquet')
GROUP BY nature
ORDER BY km DESC
LIMIT 8;┌─────────────────────┬──────────┬────────┐
│ nature │ segments │ km │
├─────────────────────┼──────────┼────────┤
│ Route à 1 chaussée │ 99343 │ 8837.0 │
│ Sentier │ 38793 │ 8167.0 │
│ Chemin │ 38202 │ 8069.0 │
│ Route empierrée │ 23118 │ 4663.0 │
│ Type autoroutier │ 727 │ 290.0 │
│ Route à 2 chaussées │ 2747 │ 281.0 │
│ Rond-point │ 5940 │ 90.0 │
│ Bretelle │ 484 │ 89.0 │
└─────────────────────┴──────────┴────────┘Trails and unpaved tracks (Sentier, Chemin, Route empierrée) add up to 20,899 km. The buildings are already planar, in EPSG:2154, so ST_Area returns m² directly, with no reprojection:
SELECT usage_1,
count(*) AS buildings,
round(avg(hauteur), 1) AS avg_height_m,
round(sum(ST_Area(geom)) / 1e6, 2) AS footprint_km2
FROM read_parquet('parquet/buildings/*.parquet')
WHERE usage_1 IS NOT NULL
GROUP BY usage_1
ORDER BY buildings DESC
LIMIT 8;┌────────────────────────┬───────────┬──────────────┬───────────────┐
│ usage_1 │ buildings │ avg_height_m │ footprint_km2 │
├────────────────────────┼───────────┼──────────────┼───────────────┤
│ Résidentiel │ 208781 │ 6.5 │ 22.92 │
│ Indifférencié │ 199543 │ 4.6 │ 11.43 │
│ Annexe │ 24221 │ 4.8 │ 1.07 │
│ Commercial et services │ 12712 │ 7.0 │ 4.88 │
│ Agricole │ 2979 │ 4.2 │ 6.43 │
│ Industriel │ 2369 │ 6.1 │ 1.69 │
│ Religieux │ 723 │ 9.5 │ 0.16 │
│ Sportif │ 318 │ 7.3 │ 0.26 │
└────────────────────────┴───────────┴──────────────┴───────────────┘Two maps
The DuckDB result converts to a GeoDataFrame in one call:
df = con.sql("SELECT importance, geom FROM read_parquet('parquet/roads/*.parquet')").df()
roads = gpd.GeoDataFrame(
df.drop(columns="geom"),
geometry=shapely.from_wkb(df["geom"].map(bytes)),
crs=4326,
).to_crs(2154)
The A9 motorway and the coastal plain around Perpignan concentrate the primary network.
For the second map, DuckDB filters the buildings to a 3.2 km square around the center of Perpignan, so that only that window crosses into Python:
SELECT hauteur, geom FROM read_parquet('parquet/buildings/*.parquet')
WHERE ST_Intersects(geom, ST_MakeEnvelope(689834, 6174971, 693034, 6178171));
The historic core, the dark cluster north of the citadel, has a median height of 12.5 m, against 7.0 m across the rest of the window.
Notes and limitations
- One SRID per column: the CRS is taken from the first non-NULL geometry and applied to the whole column. Mixed SRIDs within a column are not detected.
- Curved geometries (
CircularString,CompoundCurve,CurvePolygon) convert to the ISO WKB curve types, but reader support for those varies. BD TOPO carries only linestrings and polygons, which every reader understands.
Drawing the maps
Both maps use GeoPandas' matplotlib backend over a semi-transparent, label-free Stamen watercolor basemap (CC BY 4.0, data © OpenStreetMap contributors), pulled in with contextily. The road network is drawn from least to most important, so that the primary roads land on top, and a manual legend labels each class:
import contextily as cx
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
# Stamen retired their free tiles; this is a key-free archive mirror.
WATERCOLOR = "https://watercolormaps.collection.cooperhewitt.org/tile/watercolor/{z}/{x}/{y}.jpg"
COLORS = {
"1": "#0d366b", "2": "#184f95", "3": "#256abf",
"4": "#3987e5", "5": "#86b6ef", "6": "#c3c2b7",
}
WIDTHS = {"1": 1.8, "2": 1.3, "3": 0.9, "4": 0.55, "5": 0.3, "6": 0.2}
LABELS = {
"1": "1 - primary", "2": "2", "3": "3", "4": "4",
"5": "5 - local", "6": "paths & trails",
}
fig, ax = plt.subplots(figsize=(9, 9))
for imp in ["6", "5", "4", "3", "2", "1"]:
roads[roads.importance == imp].plot(ax=ax, color=COLORS[imp], linewidth=WIDTHS[imp], zorder=2)
cx.add_basemap(ax, crs=roads.crs, source=WATERCOLOR, alpha=0.3, attribution=False, zoom=11, zorder=1)
handles = [
Line2D([], [], color=COLORS[i], linewidth=max(WIDTHS[i], 1.2), label=LABELS[i])
for i in ["1", "2", "3", "4", "5", "6"]
]
ax.legend(
handles=handles, title="Road importance", loc="upper left",
frameon=False, labelcolor="#52514e", fontsize=9, title_fontsize=10,
)
ax.set_title("Road network of Pyrénées-Orientales - BD TOPO",
color="#0b0b0b", fontsize=14, pad=12)
ax.set_axis_off()
fig.tight_layout()
fig.savefig("img/roads_importance.png", dpi=150, bbox_inches="tight")The buildings are shaded by height on a continuous blue ramp, framed to the 3.2 km Perpignan window, over the same watercolor tiles:
from matplotlib.colors import LinearSegmentedColormap
RAMP = [
"#cde2fb", "#b7d3f6", "#9ec5f4", "#86b6ef", "#6da7ec", "#5598e7",
"#3987e5", "#2a78d6", "#256abf", "#1c5cab", "#184f95", "#104281", "#0d366b",
]
cmap = LinearSegmentedColormap.from_list("blue_ramp", RAMP)
fig, ax = plt.subplots(figsize=(9, 9))
buildings.plot(
ax=ax, column="hauteur", cmap=cmap, vmin=0, vmax=15,
missing_kwds={"color": "#e1e0d9"},
legend=True,
legend_kwds={"label": "building height (m)", "shrink": 0.55, "extend": "max"},
zorder=2,
)
ax.set_xlim(689834, 693034)
ax.set_ylim(6174971, 6178171)
cx.add_basemap(ax, crs=buildings.crs, source=WATERCOLOR, alpha=0.3, attribution=False, zoom=15, zorder=1)
ax.set_title("Buildings of central Perpignan - BD TOPO",
color="#0b0b0b", fontsize=14, pad=12)
ax.set_axis_off()
fig.tight_layout()
fig.savefig("img/perpignan_buildings.png", dpi=150, bbox_inches="tight")Package versions
Python : 3.13
SQL Server : 2022 (Docker)
FastBCP : 1.1.2
ArrowTDS : 0.5.17
geopandas : 1.1.4
shapely : 2.1.2
pyarrow : 24.0.0
duckdb : 1.5.4
matplotlib : 3.11.0
contextily : 1.7.0
