1

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!

4

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' – Saeed Soltani Moghadam Oct 30 '16 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 – Logan Byers Oct 30 '16 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' – Saeed Soltani Moghadam Oct 31 '16 at 5:15
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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.