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

Answer 2

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

def d2k(degrees, radius=6371):

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

    # Read raster data

    ds = gdal.Open(raster, 0)
    rb = ds.ReadAsArray()

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