Read and Plot a GeoTiff File Using Python

Read and Plot a GeoTiff File Using Python

Postby yiboj » Tue Aug 22, 2017 12:58 pm

GeoTiff file format is a public domain metadata standard which allows georeferencing information to be embedded within a TIFF (Tagged Image File Format) file. TIFF is a widely accepted versatile raster data format and is a standard image file format for GIS applications.

The python GDAL package is a translator library for raster and vector geospatial data formats, and it includes a variety of useful command line utilities for data translation and processing too. User can download the python package from

In this post, we use the python GDAL package to read and plot a GeoTiff file download from the LP DAAC Global Data Explorer (USGS) (

Code: Select all
import os
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from osgeo import gdal
from osgeo import osr
import numpy as np
import math

#Plot setup
fig= plt.figure(figsize=(10,8))

ax = plt.subplot(111,aspect = 'equal')
plt.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9, wspace=0, hspace=0)

#Map setup
map = Basemap(resolution='f', projection='cyl', llcrnrlon=-75, llcrnrlat=-5, urcrnrlon=-46, urcrnrlat=0)

parallels = np.arange(-50,50,1.)
meridians = np.arange(0,360,1)

map.drawparallels(parallels,labels=[1,0,0,0],color='w', fontsize=10, fontweight='bold')
meri = map.drawmeridians(meridians,labels=[0,0,0,1],color='w', fontsize=10, fontweight='bold')

#Load colormap and setup elevation contour levels
clevs = np.linspace(0, 305, 21)

# Load GeoTiff data
raster = 'geo_dem_v2_52-50.tif'

ds = gdal.Open(raster)

#Get the dimentions of column and row
nc   = ds.RasterXSize
nr   = ds.RasterYSize

#Read elevation data
data = ds.ReadAsArray()

#Get Longitude and Latitude info
geotransform = ds.GetGeoTransform()
xOrigin      = geotransform[0]
yOrigin      = geotransform[3]
pixelWidth   = geotransform[1]
pixelHeight  = geotransform[5]

#Generate Longitude and Latitude array
lons = xOrigin + np.arange(0, nc)*pixelWidth
lats = yOrigin + np.arange(0, nr)*pixelHeight

#Contour plot
x, y = map(*np.meshgrid(lons, lats))
cs=map.contourf(x, y, data, clevs, cmap=cmap)

map.drawparallels(parallels,labels=[1,0,0,0],color='k', fontsize=10, fontweight='bold')
meri = map.drawmeridians(meridians,labels=[0,0,0,1],color='k', fontsize=10, fontweight='bold')

cb = map.colorbar(cs, 'bottom', size='5%', pad='10%')

cb.set_label('Elevation (m)', fontsize=12, fontweight='bold')

Here is the image produced by the code:

Elevation Map
gdal_map.png (119.68 KiB) Viewed 10065 times
Posts: 130
Joined: Mon Mar 30, 2015 11:22 am

Return to Data Recipes