GDAL and OGR for Vectors: Inspect, Convert, Reproject, Filter, Rasterize and Polygonize

GDAL handles rasters. OGR is the part of the same library for vectors, and the two are installed and used together. In Python both come from osgeo: gdal for rasters and the tools, ogr for features and geometries. GeoPandas uses this machinery to read and write files, so what you learn here explains what happens underneath.

Why use GDAL for vectors?

GeoPandas is easier for analysis. GDAL and OGR are better for fast conversion of large files, for command line and batch jobs, and for the rasterize and polygonize steps that connect vectors and rasters.

Task	GDAL/OGR	GeoPandas
Inspect a file without loading it	`ogrinfo`	`read_file` loads it
Convert format	`ogr2ogr`	`read_file` then `to_file`
Filter by attribute	`-where`	boolean filter
Polygons to raster	`gdal_rasterize`	needs another library
Raster to polygons	`gdal_polygonize`	needs another library

1. Build sample data

Four district polygons with an attribute, written to a GeoPackage with GeoPandas (Part 10). The coordinates are in EPSG:32643 and the data is invented.

# Step 1: import the tools
import geopandas as gpd
import numpy as np
from pathlib import Path
from osgeo import gdal, ogr, osr          # GDAL: gdal = rasters, ogr = vectors, osr = coordinate systems
from shapely.geometry import box

# Step 2: ask GDAL to stop with a clear Python error when something fails
# (without this, it can fail quietly and return nothing)
gdal.UseExceptions()
ogr.UseExceptions()

# Step 3: make a folder called "vectors" for all the files in this part
Path("vectors").mkdir(exist_ok=True)

# Step 4: make four district rectangles. We use GeoPandas here only to create test data.
X0, Y0 = 700000, 1400000                  # a corner point, in metres
districts = gpd.GeoDataFrame(
    {
        "name": ["North", "South", "East", "West"],
        "zone_code": [1, 2, 3, 2],                     # a number for each zone
        "density": [320.5, 85.0, 140.2, 60.7],         # people per km2
    },
    geometry=[                                         # box(left, bottom, right, top)
        box(X0, Y0 + 30000, X0 + 60000, Y0 + 60000),
        box(X0, Y0, X0 + 40000, Y0 + 30000),
        box(X0 + 40000, Y0, X0 + 60000, Y0 + 30000),
        box(X0 - 30000, Y0, X0, Y0 + 60000),
    ],
    crs="EPSG:32643",
)

# Step 5: save it as a GeoPackage file that GDAL will read in the next steps
districts.to_file("vectors/districts.gpkg", layer="districts", driver="GPKG")
print(len(districts), "districts saved")

4 districts saved

2. Inspect with OGR

This is the equivalent of ogrinfo.

# Step 1: open the file. A file (a "datasource") can hold several layers.
ds = ogr.Open("vectors/districts.gpkg")
print("layers:", ds.GetLayerCount())

# Step 2: take the first layer (counting starts at 0) and ask it basic questions
layer = ds.GetLayer(0)
print("name:", layer.GetName())
print("features:", layer.GetFeatureCount())                       # how many rows
print("geometry:", ogr.GeometryTypeToName(layer.GetGeomType()))   # polygon, point, line...

# Step 3: the coordinate system, as an authority and a code, e.g. EPSG 32643
srs = layer.GetSpatialRef()
print("crs:", srs.GetAuthorityName(None), srs.GetAuthorityCode(None))
print("extent:", layer.GetExtent())                               # min x, max x, min y, max y

# Step 4: list the attribute columns (fields). The "definition" describes the table structure.
definition = layer.GetLayerDefn()
for i in range(definition.GetFieldCount()):                       # repeat once for each field
    field = definition.GetFieldDefn(i)
    print(" field:", field.GetName(), field.GetTypeName())

layers: 1
name: districts
features: 4
geometry: Polygon
crs: EPSG 32643
extent: (670000.0, 760000.0, 1400000.0, 1460000.0)
 field: name String
 field: zone_code Integer64
 field: density Real

A data source (ds) holds one or more layers. A GeoPackage can hold many, a Shapefile one. Each layer has a geometry type, a CRS, an extent given as (minx, maxx, miny, maxy) and a list of fields. Note the extent order of OGR, which differs from GeoPandas' total_bounds. Loop over the features to see their attributes and geometry:

# A layer can be looped over: each pass gives one feature (one row with its shape)
for feature in layer:
    geom = feature.GetGeometryRef()                               # the shape of this feature
    # GetArea is in square metres (the CRS is in metres), so divide by 1,000,000 for km2
    print(feature.GetFID(), feature["name"], feature["density"], round(geom.GetArea() / 1e6), "km2")

1 North 320.5 1800 km2
2 South 85.0 1200 km2
3 East 140.2 600 km2
4 West 60.7 1800 km2

The command line version of all of this is ogrinfo vectors/districts.gpkg districts -so, where -so means summary only. Without -so it prints every feature.

3. Convert formats

gdal.VectorTranslate is ogr2ogr: source in, destination out, with options.

# VectorTranslate converts a vector file into another format (like ogr2ogr on the command line).
# The pattern is: (output file, input file, format=...)
gdal.VectorTranslate("vectors/districts.geojson", "vectors/districts.gpkg", format="GeoJSON")
gdal.VectorTranslate("vectors/districts.shp", "vectors/districts.gpkg", format="ESRI Shapefile")

# Open each new file and check that it kept the rows, shape type and field names
for name in ["districts.geojson", "districts.shp"]:
    data = ogr.Open(f"vectors/{name}")
    lyr = data.GetLayer(0)
    field_names = [lyr.GetLayerDefn().GetFieldDefn(i).GetName()
                   for i in range(lyr.GetLayerDefn().GetFieldCount())]
    print(name, lyr.GetFeatureCount(), ogr.GeometryTypeToName(lyr.GetGeomType()), field_names)

districts.geojson 4 Polygon ['name', 'zone_code', 'density']
districts.shp 4 Polygon ['name', 'zone_code', 'density']

After converting, verify the feature count, geometry type, fields and CRS. Check the CRS of the result and do not assume it:

# Which coordinate system did the GeoJSON file keep? GDAL does not reproject on its own,
# so the metre-based code of the original file is carried over.
data = ogr.Open("vectors/districts.geojson")
print(data.GetLayer(0).GetSpatialRef().GetAuthorityCode(None))

GDAL writes the CRS into the GeoJSON file and keeps EPSG:32643. The GeoJSON standard (RFC 7946) expects EPSG:4326, and many web tools assume it. To follow the standard, reproject while converting, as in the next step, or add the option layerCreationOptions=["RFC7946=YES"], which writes coordinates in EPSG:4326.

The command line forms are ogr2ogr -f GeoJSON districts.geojson districts.gpkg and ogr2ogr -f "ESRI Shapefile" districts.shp districts.gpkg. Shapefile field names are limited to 10 characters, and it is made of several files (.shp, .shx, .dbf, .prj) that must stay together.

4. Reproject

Add dstSRS to move to another CRS. This is ogr2ogr -t_srs.

# Step 1: convert the coordinates while copying. dstSRS = "destination spatial reference system".
gdal.VectorTranslate(
    "vectors/districts_wgs84.gpkg", "vectors/districts.gpkg",
    format="GPKG", dstSRS="EPSG:4326", layerName="districts",
)

# Step 2: open the result and confirm the code and the new extent (now in degrees)
data = ogr.Open("vectors/districts_wgs84.gpkg")
lyr = data.GetLayer(0)
print(lyr.GetSpatialRef().GetAuthorityCode(None))
print([round(v, 4) for v in lyr.GetExtent()])

# Step 3: do the same with GeoPandas and compare. The numbers differ only in order:
# GDAL gives min x, max x, min y, max y; GeoPandas gives min x, min y, max x, max y.
via_geopandas = gpd.read_file("vectors/districts.gpkg").to_crs(4326)
print(via_geopandas.total_bounds.round(4).tolist())

4326
[76.5654, 77.3989, 12.6534, 13.2019]
[76.5654, 12.6534, 77.3989, 13.2019]

The two results agree: GDAL gives (minx, maxx, miny, maxy) and GeoPandas gives (minx, miny, maxx, maxy), the same four numbers in a different order. Both tools use the same underlying coordinate engine (PROJ). Use whichever fits the rest of your script.

5. Filter features

A where clause is a SQL condition on the attributes. It is the ogr2ogr -where option.

# Step 1: copy only the rows where density is above 100. The "where" text works like SQL.
gdal.VectorTranslate(
    "vectors/dense.gpkg", "vectors/districts.gpkg",
    format="GPKG", where="density > 100", layerName="dense",
)

# Step 2: open the result and list the names that survived
data = ogr.Open("vectors/dense.gpkg")
lyr = data.GetLayer(0)
print(lyr.GetFeatureCount(), [f["name"] for f in lyr])

# Step 3: the same filter in GeoPandas, to see that both give the same answer
same = gpd.read_file("vectors/districts.gpkg")
print(same[same["density"] > 100]["name"].tolist())

2 ['North', 'East']
['North', 'East']

The GeoPandas filter in the last two lines is the Part 7 pattern and gives the same answer. For large files GDAL does the filtering while reading, so it never loads the rest. Use SQLStatement="SELECT name, density FROM districts WHERE zone_code = 2", SQLDialect="SQLITE" for selections that need SQL, such as choosing columns.

6. Rasterize: polygons to cells

gdal.Rasterize burns polygon values into a raster grid. You choose the grid.

# Rasterize turns polygons into a grid of cells (a raster).
# Every cell takes the "density" value of the polygon it falls in.
gdal.Rasterize(
    "vectors/density.tif", "vectors/districts.gpkg",   # output raster, input vector
    attribute="density",                               # which column becomes the cell value
    outputBounds=[670000, 1390000, 770000, 1470000],   # the area to cover: min x, min y, max x, max y
    xRes=1000, yRes=1000,                              # each cell is 1000 m x 1000 m
    outputType=gdal.GDT_Float64,                       # cells hold decimal numbers
    noData=-9999,                                      # cells outside every polygon get this value
    outputSRS="EPSG:32643",
)

# Read the cell values into a NumPy array, which is Python's table of numbers
ds = gdal.Open("vectors/density.tif")
arr = ds.GetRasterBand(1).ReadAsArray()                # band 1 = the first (and only) layer of values
print(arr.shape)                                       # (rows, columns)
print([round(float(v), 1) for v in np.unique(arr)])    # the distinct values found

(80, 100)
[-9999.0, 60.7, 85.0, 140.2, 320.5]

outputBounds is the extent as minx, miny, maxx, maxy, and xRes, yRes are the cell size in CRS units, here 1 km. Every cell takes the density of the polygon that covers its centre, and cells outside every polygon get the NoData value. The bounds are a little larger than the districts, so the raster has 80 rows by 100 columns and the unique values are the four densities plus the NoData value for the empty margin. If polygons overlap, the later one in the file overwrites the earlier one. The choice of cell size matters: a coarser grid makes the borders more jagged. The command line form is gdal_rasterize -a density -tr 1000 1000 -te 670000 1390000 770000 1470000 -a_nodata -9999 districts.gpkg density.tif.

7. Polygonize: cells to polygons

gdal.Polygonize does the reverse. It groups neighbouring cells with the same value into polygons. Make a class raster first, and then convert it back.

# Polygonize goes the other way: grid cells with the same value become polygons.

# Step 1: make a small raster from the zone_code column (whole numbers, 0 = no data)
gdal.Rasterize(
    "vectors/zones.tif", "vectors/districts.gpkg",
    attribute="zone_code",
    outputBounds=[670000, 1390000, 770000, 1470000],
    xRes=1000, yRes=1000,
    outputType=gdal.GDT_Int16,
    noData=0,
    outputSRS="EPSG:32643",
)

# Step 2: open it and take its first band
src = gdal.Open("vectors/zones.tif")
band = src.GetRasterBand(1)

# Step 3: prepare an empty output file with one layer and one integer field for the zone number
drv = ogr.GetDriverByName("GPKG")                      # the driver knows how to write GeoPackages
out = drv.CreateDataSource("vectors/zones_poly.gpkg")
srs = osr.SpatialReference()
srs.ImportFromWkt(src.GetProjection())                 # reuse the raster's coordinate system
out_layer = out.CreateLayer("zones", srs, ogr.wkbPolygon)
out_layer.CreateField(ogr.FieldDefn("zone_code", ogr.OFTInteger))

# Step 4: do the conversion
# Arguments: input band, mask band, output layer, index of the field that stores the value (0 = first)
gdal.Polygonize(band, band, out_layer, 0)      # the second band is a mask: NoData cells are skipped
out = None                                     # setting it to None closes and saves the file

# Step 5: read the result back with GeoPandas and measure the polygons
polys = gpd.read_file("vectors/zones_poly.gpkg")
polys["area_km2"] = (polys.geometry.area / 1e6).round(0)
print(polys.sort_values("zone_code")[["zone_code", "area_km2"]])
print(len(polys), "polygons")

   zone_code  area_km2
0          1    1800.0
1          2    3000.0
2          3     600.0
3 polygons

Each polygon gets the cell value in the field you point to, and the last number in the call is the index of that field (0 is the first one). Using the band as its own mask skips the NoData cells, which have the value 0. Four districts become three polygons, because South and West touch and share zone code 2, so Polygonize merges them and the border between them disappears. Polygonize groups touching cells of equal value, whatever the original features were. The polygon edges follow the cell edges, so they are stair-stepped. The command line tool is gdal_polygonize.py zones.tif -f GPKG zones_poly.gpkg zones zone_code.

8. Automate a batch

A loop with a log turns the tools into a repeatable job. This converts every GeoPackage in a folder to GeoJSON and writes a record of what it did.

# A mini batch job: convert every GeoPackage in the folder to GeoJSON and keep a log of what happened.
import csv

# Step 1: create an output folder and an empty list to collect one log row per file
out_dir = Path("vectors/geojson_out")
out_dir.mkdir(exist_ok=True)
log = []

# Step 2: go through every .gpkg file, in alphabetical order
for path in sorted(Path("vectors").glob("*.gpkg")):
    target = out_dir / (path.stem + ".geojson")        # same name, new extension
    try:                                               # try the conversion...
        gdal.VectorTranslate(str(target), str(path), format="GeoJSON", dstSRS="EPSG:4326")
        data = ogr.Open(str(target))
        count = data.GetLayer(0).GetFeatureCount()
        data = None
        log.append({"input": path.name, "output": target.name, "features": count, "status": "ok"})
    except RuntimeError as error:                      # ...and if it fails, note the reason and carry on
        log.append({"input": path.name, "output": "", "features": 0, "status": str(error)[:60]})

# Step 3: write the log to a CSV file (a table you can open in Excel)
with open(out_dir / "log.csv", "w", newline="") as f:
    writer = csv.DictWriter(f, fieldnames=["input", "output", "features", "status"])
    writer.writeheader()                               # the column names
    writer.writerows(log)                              # one row per file

# Step 4: print the same log on screen
for row in log:
    print(row["input"], "->", row["output"], row["features"], row["status"])

dense.gpkg -> dense.geojson 2 ok
districts.gpkg -> districts.geojson 4 ok
districts_wgs84.gpkg -> districts_wgs84.geojson 4 ok
zones_poly.gpkg -> zones_poly.geojson 3 ok

Each file is converted, counted and logged, and a failure is recorded without stopping the loop. The same pattern works for rasters. The loop and the log follow the same pattern as the file inventory in Part 5 and the city cleaner in Part 6.

Run this part as a notebook

The code on this page is also a Jupyter notebook. It has the theory in text cells, then the code, then guided practice, an independent assignment and hidden solutions, so you can try each task before you look. Download the Part 13 notebook, open it in JupyterLab and run the cells from top to bottom. Part 1 shows how to install and start Jupyter. This part needs the GDAL Python bindings, which come with QGIS. The notebook was run on QGIS's own Python, and if GDAL is installed in another way, check that from osgeo import gdal works before you start.

Practice

Run ogrinfo -so on a vector file of your own and list its layers, CRS and fields.
Convert a Shapefile to GeoPackage and compare the field names.
Filter a layer with where using two conditions joined by AND.
Rasterize a polygon layer at two cell sizes and compare the raster sizes.
Polygonize a land cover raster and use GeoPandas to total the area per class.

Common mistakes

Reading the OGR extent as (minx, miny, maxx, maxy). The order is (minx, maxx, miny, maxy).
Not closing a data source you wrote (out = None), so the file is incomplete.
Chaining ogr.Open(path).GetLayer(0) in one line. The data source is discarded and the layer stops working, with a TypeError. Keep the data source in a variable for as long as you use its layer.
Mixing up the CRS of the vector and the raster when rasterizing.
Choosing a cell size that is much coarser than the smallest polygon.
Assuming a Shapefile keeps long field names or Unicode.
Polygonizing a continuous raster, which creates one polygon per cell.

Key takeaways

Next: SQL basics for GIS data.