[GDAL-Python] Create & save geotiff using gdal in python

In this tutorial, we will see how to create & save an numpy array data as geotiff.

1. Import the required packages — gdal, matplotlib & numpy

from osgeo import gdal, osr
import matplotlib.pyplot as plt
import numpy as np

2. Let’s create a sample raster data to save it as a geotiff. We have generated an image data with 2 bands and is visualized below.

#sample image data
data_arr = np.array([[
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,5,4,0,0,0,5,4,0,0,4,4,3,2,0],
    [0,4,0,4,0,4,0,4,0,5,0,0,0,0,0],
    [0,3,0,0,3,0,0,3,0,4,0,0,0,0,0],
    [0,2,0,0,0,0,0,2,0,4,0,0,0,0,0],
    [0,1,0,0,0,0,0,1,0,0,3,2,1,1,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
],
[
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,1,1,0,0,0,1,1,0,0,1,1,1,1,0],
    [0,1,0,1,0,1,0,1,0,1,0,0,0,0,0],
    [0,1,0,0,1,0,0,1,0,1,0,0,0,0,0],
    [0,1,0,0,0,0,0,1,0,1,0,0,0,0,0],
    [0,1,0,0,0,0,0,1,0,0,1,1,1,1,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
]])
#visualising sample data
plt.figure(figsize = (10, 10))
plt.imshow(data_arr[0], cmap='gray')
plt.show()
plt.figure(figsize = (10, 10))
plt.imshow(data_arr[1], cmap='gray')

3. Now, we have to define few parameters — Image geographical extent, resolution, width & height.

# width & height for any band
no_of_bands = data_arr.shape[0]
height = data_arr[0].shape[0]
width = data_arr[0].shape[1]


#sample geoextent & resolution - Update it according to your requirement
x1 = 139.7524088
y1 = 35.9806627
x2 = 139.7557615
y2 = 35.9790435
x_res = (x1 - x2)/width
y_res = (y1 - y2)/width
print(no_of_bands, x1, y1, x2, y2, x_res, y_res, width, height)

4. Now that our sample data is ready, let’s start creating a tiff image. First, load the GTiff driver from gdal.

driver = gdal.GetDriverByName('GTiff')

5. Create the target dataset.

out_ds = driver.Create('./output.tif', width, height, no_of_bands, gdal.GDT_Float32)

6. Update the GeoTransform of the dataset.

out_ds.SetGeoTransform((x1, x_res, 0, y1, 0, y_res))

7. Write data array into dataset bands.

for band_no in range(0, no_of_bands):
    band = out_ds.GetRasterBand(band_no + 1)
    band.WriteArray(data_arr[band_no])
    band.FlushCache()

8. Update the spatial reference system of the dataset & close the dataset.

out_srs = osr.SpatialReference()
out_srs.ImportFromEPSG(4326)
out_ds.SetProjection(out_srs.ExportToWkt())
out_ds = None

9. Read & Visualise the data using rasterio.

# read and visualize the saved geotiff
import rasterio
from rasterio.plot import show
src = rasterio.open('./output.tif')
fig, ax = plt.subplots(1, figsize=(12, 12))
show(src, ax=ax)