Geospatial · Practical guide
Read & write GeoTIFF
Move between georeferenced raster files and NumPy arrays while keeping the spatial metadata intact.
Start learning01GeoTIFF02NumPy array03GeoTIFF
What you’ll learn
- Read raster values and spatial metadata.
- Write multiband arrays back to a GeoTIFF.
- Check array shape, coordinate system, and geotransform.
Before you begin Python basics
- A Python environment with NumPy and GDAL.
- A local GeoTIFF for the usage example.
- Arrays use rows × columns × bands; a single band needs a band axis.
01
Read & write helpers
The reader returns the array, geotransform, projection, and raster dimensions. The writer expects a three-dimensional array and writes Float32 bands.
- Array shape
- Multiband reads are rearranged from bands × rows × columns to rows × columns × bands.
- Spatial metadata
- Pass the original geotransform and projection to the writer to retain the raster location.
import numpy as np
from osgeo import gdal, ogr, gdalconst
def read_tif(tif_file):
"""
Parameters:
------------------------------
- tif_file: the full path of the GeoTIFF data (e.g., "/data/imagery.tif").
- im_data: (rows, cols, bands)
"""
dataset = gdal.Open(tif_file)
cols = dataset.RasterXSize
rows = dataset.RasterYSize
im_proj = (dataset.GetProjection())
im_Geotrans = (dataset.GetGeoTransform())
im_data = dataset.ReadAsArray(0, 0, cols, rows)
if im_data.ndim == 3:
im_data = np.moveaxis(dataset.ReadAsArray(0, 0, cols, rows), 0, -1)
dataset = None
return im_data, im_Geotrans, im_proj, rows, cols
def array_to_geotiff(array, output_path, geo_transform, projection, band_names=None):
"""
Parameters:
------------------------------
- array: the image array that need to be saved (rows, cols, bands).
- output_path: the full path that save the GeoTIFF data (e,g., "/data/saved_imagery.tif").
- geo_transform: geotransform of data.
- projection: projection of data.
- band_names: the band names as list.
"""
rows, cols, num_bands = array.shape
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create(output_path, cols, rows, num_bands, gdal.GDT_Float32)
dataset.SetGeoTransform(geo_transform)
dataset.SetProjection(projection)
for band_num in range(num_bands):
band = dataset.GetRasterBand(band_num + 1)
band.WriteArray(array[:, :, band_num])
band.FlushCache()
if band_names:
band.SetDescription(band_names[band_num])
dataset = None
band = None
return02
Use the helpers
Save the helpers above as read_write_geotiff.py, place a raster in your data folder, and replace the input and output filenames.
from read_write_geotiff import read_tif, array_to_geotiff
array, transform, projection, rows, cols = read_tif("data/imagery.tif")
if array.ndim == 2:
array = array[..., None]
array_to_geotiff(array, "data/imagery_copy.tif", transform, projection)
print(rows, cols, array.shape)
03
Check the output
After writing, inspect the output alongside the input. A matching array shape alone does not confirm that the coordinates were retained.
- Confirm rows, columns, and band count.
- Compare the projection and six geotransform values.
- Check data type and NoData behavior for your analysis.
- Reopen the output and compare raster values.