# Plotting cross section using python, gdal, tiff?

How to plot a cross section between two pairs of lat/long using a tiff file in python?

At first glance, it seems i have to calculate points position between these two desired points and then looking for nearest points in tiff to get their elevation.

Until now i could:

ds = gdal.Open('cliped.tif')

but the problem is data point has no projection so i can not search for lat/long pairs!

You have the start of the solution to your problem. Here is the basic workflow:

1. Open raster, get raster coordinate system
2. Get points as x,y coordinates in same coordinate system
3. Get a function that transforms x,y into raster i,j coordinates
4. Extract values along the line segment between original points
5. Plot

The code follows:

import gdal
import osr
import matplotlib.pyplot as plt

raster = 'input.tif'
ds = gdal.Open(raster, 0)
rb = ds.GetRasterBand(1)

gt = ds.GetGeoTransform() # maps i,j to x,y
# try-except block to handle different output of InvGeoTransform with gdal versions
try:
inv_gt_success, inverse_gt = gdal.InvGeoTransform(gt) # maps x,y to i,j
except:
inverse_gt = gdal.InvGeoTransform(gt) # maps x,y to i,j

sr_ds = osr.SpatialReference()  # spatial reference of the dataset
sr_ds.ImportFromWkt(ds.GetProjection())
sr_wgs84 = osr.SpatialReference()  # spatial reference of WGS84
sr_wgs84.SetWellKnownGeogCS('WGS84')
ct = osr.CoordinateTransformation(sr_wgs84, sr_ds)

pt0 = (30., 60.) # lat lon
pt1 = (35., 65.) # lat lon
pt0_ds = ct.TransformPoint(*pt0)  # x,y
pt1_ds = ct.TransformPoint(*pt1)  # x,y

num_pts = 100
dx = (pt1_ds[0] - pt0_ds[0]) / num_pts
dy = (pt1_ds[1] - pt0_ds[1]) / num_pts

raster_vals = []  # stores the extracted values

for i in range(num_pts):
point = (pt0_ds[0] + i * dx, pt0_ds[1] + i * dy)

pix_x = int(inverse_gt[0] + inverse_gt[1] * point[0] +
inverse_gt[2] * point[1])
pix_y = int(inverse_gt[3] + inverse_gt[4] * point[0] +
inverse_gt[5] * point[1])
val = rb.ReadAsArray(pix_x, pix_y, 1, 1)[0,0]
raster_vals.append(val)

plt.plot(raster_vals)
plt.show()
• Many Many Thanks Dear Logan ;), But it seems something wrong when i run your code. it gives the following error: Traceback (most recent call last): File "/home/saeed/Program/seismo/WOR/PhD/1D_velocity_model/PyHYP/t.py", line 34, in <module> pix_x = int(inverse_gt[0] + inverse_gt[1] * point[0] + TypeError: can't multiply sequence by non-int of type 'float' Oct 30, 2016 at 17:43
• @SaeedSoltaniMoghadam The example is updated. Different versions of gdal have different outputs for InvGeoTransform(). Your version outputs two vairables, a success flag and the inverse geo transform Oct 30, 2016 at 18:25
• I did it but i got no answer. Still give the following errors! Traceback (most recent call last): File "/home/saeed/Program/seismo/WOR/PhD/1D_velocity_model/PyHYP/cross.py", line 40, in <module> val = rb.ReadAsArray(pix_x, pix_y, 1, 1)[0,0] TypeError: 'NoneType' object has no attribute 'getitem' Oct 31, 2016 at 5:15

Thanks to @Logan, here is my own code and now it works:

import gdal
import osr
import matplotlib.pyplot as plt
from numpy import arange, array, nonzero, in1d, abs, linspace, sqrt, pi

# Conver Degree to Kilometer (Thanks to Obspy!)

return degrees * (2.0 * radius * pi / 360.0)

# Find Nearest point in array

def find_nearest(array,value):

idx = (abs(array-value)).argmin()

return idx

def plot_cross_elv(raster=None, A=None, B=None, np=1000):

ds = gdal.Open(raster, 0)

# Get start/end point of the map in lon/lat

coords = ds.GetGeoTransform()

# Get lon/lat and elevation

lon = array([coords[0]+coords[1]*i for i in range(rb.shape[1])])
lat = array([coords[3]+coords[5]*i for i in range(rb.shape[0])])
elv = rb.flatten()

# Make pairs of (lon,lat)

points = array([array([x,y]) for y in lat for x in lon])

# Define Profile using lon/lat

ABx   = linspace(A[0],B[0],np)
ABy   = linspace(A[1],B[1],np)

profile = array([array([x,y]) for x,y in zip(ABx, ABy)])

# Find Nearest points of profile in Map and Prepare for plot

cross_dis = []
cross_elv = []

for p in profile:

lon_ind = find_nearest(lon,p[0])
lat_ind = find_nearest(lat,p[1])

x = lon[lon_ind]
y = lat[lat_ind]
d = sqrt(d2k(A[0]-x)**2+d2k(A[1]-y)**2)
h = rb[lat_ind][lon_ind]

cross_dis.append(d)
cross_elv.append(h)

ax = plt.subplot(111)
ax.plot(cross_dis,cross_elv,linewidth=2,color='k')
ax.set_xlabel('Distance [km]')
ax.set_ylabel('Elevation [m]')
ax.fill_between(cross_dis, 0, cross_elv, color='bisque')
ax.set_xlim(min(cross_dis), max(cross_dis))
ax.locator_params(axis='x',nbins=6)
ax.locator_params(axis='y',nbins=6)

plt.show()

# RUN

plot_cross_elv(raster='cliped.tif', A=(50.0, 35.0), B=(52.0, 36.5), np=1000)