GDAL (Geospatial Data Abstraction Library)

Warnings

• gotcha The `pip install GDAL` command only installs the Python bindings. The core GDAL C/C++ library and its dependencies (like PROJ and GEOS) must be installed separately on your operating system. Failure to do so will result in build errors during pip installation or a non-functional Python package. Consult GDAL's official documentation or OS-specific package managers (e.g., `apt`, `yum`, `brew`, `conda`) for proper system-level installation.
  Fix: Ensure the GDAL C/C++ library is installed on your system *before* attempting `pip install GDAL`.
• breaking GDAL 3.0 introduced significant changes, especially regarding Coordinate Reference Systems (CRS) and axis order. It standardized on 'east/north' (longitude/latitude) axis order for geographic CRS, aligning with ISO 19111, whereas GDAL 2.x often used 'north/east'. This can lead to incorrect coordinate transformations or interpretations if not explicitly managed, especially when working with older datasets or applications.
  Fix: Be mindful of axis order. For `osr.SpatialReference` objects, use `SetAxisMappingStrategy()` to explicitly control how axis order is handled (e.g., `osr.OAMS_TRADITIONAL_GIS_ORDER` for Lat/Lon or `osr.OAMS_AUTHORITY_COMPLIANT` for Lon/Lat). Test existing code paths that perform coordinate transformations.
• gotcha GDAL datasets consume system resources and often hold open file handles. It is crucial to explicitly close datasets when you are finished with them to prevent resource leaks and ensure data integrity (e.g., for writing operations).
  Fix: Set the dataset object to `None` (e.g., `dataset = None`) or allow it to go out of scope. For multiple datasets, ensure all are explicitly closed. Using `try...finally` blocks or context managers (if available for specific operations) can help manage this.
• deprecated The default output of Well-Known Text (WKT) representations for spatial references changed from WKT1 to WKT2 in GDAL 3.0. While GDAL can still parse WKT1, generating WKT1 by default is no longer the standard. This can affect interoperability with older GIS software that may not fully support WKT2.
  Fix: Use `osr.SpatialReference.ExportToWkt(options=['VERSION=WKT1_GDAL'])` to explicitly generate WKT1 if backward compatibility is required. Otherwise, ensure consuming applications can handle WKT2.
• gotcha GDAL's C++ core often uses C-style error reporting (e.g., logging to stderr, returning NULL or specific error codes) rather than throwing exceptions. While Python bindings translate some errors to exceptions, others might result in `None` being returned or a silent failure with logged warnings. This can make debugging difficult.
  Fix: Call `gdal.UseExceptions()` at the start of your script to make GDAL raise Python exceptions for most errors, allowing for more robust error handling with `try...except` blocks. Always check return values for `None` from GDAL functions.

Install

• pip install GDAL Install Python Bindings (requires system GDAL)

Imports

• gdal

  from osgeo import gdal

• ogr

  from osgeo import ogr

  For vector data operations.
• osr

  from osgeo import osr

  For spatial reference system operations.
• gdalconst

  from osgeo import gdalconst

  For GDAL constants like data types, access modes.

Quickstart

This quickstart demonstrates basic raster and vector data handling using in-memory datasets with `osgeo.gdal` and `osgeo.ogr`. It includes setting up error handling to raise Python exceptions and ensures datasets are properly released.

from osgeo import gdal, ogr, osr

# Configure GDAL to raise Python exceptions for errors
gdal.UseExceptions()

# --- Raster Quickstart (using an in-memory dataset) ---
print("--- Raster Quickstart ---")
driver = gdal.GetDriverByName("MEM")
# Create a 10x10, 1-band, Byte-type in-memory raster dataset
ds_raster = driver.Create("dummy_raster", 10, 10, 1, gdal.GDT_Byte)
ds_raster.SetProjection("EPSG:4326") # Set a dummy geographic projection
raster_band = ds_raster.GetRasterBand(1)
raster_band.WriteArray([[i * j for j in range(10)] for i in range(10)])

print(f"Raster Driver: {ds_raster.GetDriver().LongName}")
print(f"Raster Size: {ds_raster.RasterXSize}x{ds_raster.RasterYSize}")
print(f"Raster Projection: {ds_raster.GetProjectionRef()[:20]}...") # Abbreviate for brevity
print(f"First pixel value: {raster_band.ReadAsArray()[0, 0]}")

ds_raster = None # Release the dataset

# --- Vector Quickstart (using an in-memory dataset) ---
print("\n--- Vector Quickstart ---")
driver = ogr.GetDriverByName("Memory")
# Create an in-memory vector datasource
ds_vector = driver.CreateDataSource("dummy_vector")
# Create a layer
srs = osr.SpatialReference()
srs.SetFromUserInput("EPSG:4326")
layer = ds_vector.CreateLayer("points", srs, ogr.wkbPoint)

# Add a field
field_defn = ogr.FieldDefn("name", ogr.OFTString)
layer.CreateField(field_defn)

# Create a feature
feature = ogr.Feature(layer.GetLayerDefn())
feature.SetField("name", "Test Point")
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10.0, 20.0) # Lon, Lat
feature.SetGeometry(point)
layer.CreateFeature(feature)

print(f"Vector Layer Name: {layer.GetName()}")
print(f"Number of Features: {layer.GetFeatureCount()}")
# Access the first feature
first_feature = layer.GetFeature(0)
print(f"First Feature Name: {first_feature.GetField("name")}")
print(f"First Feature Geometry: {first_feature.GetGeometryRef().ExportToWkt()}")

ds_vector = None # Release the datasource