Geographical data visualization


Author: GC

Date: Aug. 13, 2016

Previous Poincaré-Birkhoff points

Next Saddle node bifurcation at infinity


The next program uses matplotlib to visualize some geographical elevation data. The data can be collected on the website http://srtm.csi.cgiar.org/ , and converted to an exploitable format through the following bash command:

gdal_translate -of XYZ srtm_05_08.tif srtm_05_08.xyz



Geographical data visualization




Source code

#! /usr/bin/env python3

"""
Programm to visualize geographical data collected from http://srtm.csi.cgiar.org/
Current example: srtm_05_08.xyz -> Hawaii islands
"""

# import matplotlib and numpy
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.mlab import griddata
from matplotlib import cm

# data extraction
file = open('/path/srtm_05_08/srtm_05_08.xyz')
data = open('/path/srtm_05_08/hawaii.txt', 'w')

for line in file:
    Ligne = file.readline()
    Ligne = Ligne.split()
    
    if len(Ligne)==3:
        x = float(Ligne[0])
        y = float(Ligne[1])
        z = float(Ligne[2])
        
        # select the data corresponding to some region
        latitude = 20.896
        longitude = -156.585
        eps = 0.1
        
        if abs(x-longitude)<eps and abs(y-latitude)<eps:
            data.write(str(x)+' '+str(y)+' '+str(z)+'\n')
            
file.close()
data.close()
print('Read data: OK')

# read data file
data = np.genfromtxt('/path/srtm_05_08/hawaii.txt')
print('Extract data: OK')

x = data[1:,0]
y = data[1:,1]
z = data[1:,2]

# correction
for i in range(len(z)):
    if z[i]<=0:z[i]=0

# grid generation
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)

# figure
fig = plt.figure()
ax = fig.gca(projection='3d')
plt.setp(ax.get_xticklabels(), visible=False)
plt.setp(ax.get_xticklines(), visible=False)
plt.setp(ax.get_xgridlines(), visible=False)
plt.setp(ax.get_yticklabels(), visible=False)
plt.setp(ax.get_yticklines(),visible=False)
plt.setp(ax.get_zticklabels(), visible=False)
plt.setp(ax.get_zticklines(),visible=False)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_zticklabels([])
ax.grid(False)
cb = cm.terrain
Z = griddata(x, y, z, xi, yi, interp='linear')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cb, vmin=0, vmax=max(z), linewidth=0.1, antialiased=True)
ax.set_zlim(0, 2*max(z))
plt.show()


Previous Poincaré-Birkhoff points

Next Saddle node bifurcation at infinity