[GIS] Plotting cross section using python, gdal, tiff

Question

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:

# Read tiff data
ds = gdal.Open('cliped.tif')
# Read it as Array
a=ds.ReadAsArray()

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

Answer 1

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()