Data analysis and visualization of digital elevation of Bangladesh

The digital elevation model is basically a 3D representation of the earth's surface or any terrain’s surface, it is well known as DEM which is mainly created from elevation data. A terrain is mathematically modeled as a function z = f(x, y) which maps each point (x, y) in a planar domain D into an elevation value f(x, y). In this view, the terrain is the graph of function f over D.DEM can be represented as a raster file [grid of squares] with each pixel value corresponding to a height above a datum or vector-based triangular irregular network (TIN).

We can say, it is an integral part of GIS and used for 1) hydrological modeling including flood simulation, delineation, and analysis of watershed and drainage network.

2) Soil erosion and sediment transport modeling, 3) delineation and study of physiographic units

4) soil and ecological studies,5) geomorphological evolution of landforms,6) and civil engineering and military applications such as site and route selection, viewshed analysis, landslide hazard assessment, and so on.

Today, I am going to show how to visualize digital elevation data, download, extraction of latitude, longitude, and height, automatic plotting from elevation data, automatic CSV file generation, making slope,

Data source: SRTM data, spatial resolution approximately 90 meters on the line of the equator, for the entire globe are available:

Data description: Product: SRTM 90m DEM Version 4
Data File Name: srtm_55_08.zip
Mask File Name: srtm_mk_55_08.zip
Latitude Min: 20 N Max: 25 N
Longitude Min: 90 E Max: 95 E
Center Point Lat: 22.5 N Long: 92.5 E

you can download it from this link: http://srtm.csi.cgiar.org/

In this analysis and visualization, I will use some python and R libraries

1. Rasterio, raster,plot3d

2. geopandas

3. pandas

4. numpy

5. matplotlib

I will use the most widely used open-access global DEM (90-meter SRTM) is almost two decades old. http://srtm.csi.cgiar.org/

Figure 1 typical Digital Elevation Representation

1: I want to plot DEM in the machine without using commercial software or GUI based software.

code :

rasterVis::levelplot(r1,
margin = list(x = FALSE,
y = TRUE),
main=’Digital Elevation of Chittagong’,
col.regions = terrain.colors(16),
xlab = list(label = “”,
vjust = -0.25),
sub = list(
label = “elevation”,
font = 1,
cex = .9,
hjust = 1.5))

this is a typical visualization of Dem with R

Now, Get slope and aspect data and plot from this DEM data.

CODE:

x <- terrain(r1, opt = c(“slope”, “aspect”), unit = “degrees”)
plot(x, main=’Slope and aspect of Chittagong’)

The slope of this DEM

the aspect of this DEM

Open the file and get details

file<- ‘C:/Users/hafez/personal/GIS DATA/srtm_55_08/srtm_55_08/srtm_55_08.tif’
r1= raster(file)

#get detials

r1, class : RasterLayer, dimensions : 6000, 6000, 3.6e+07 (nrow, ncol, ncell), resolution : 0.0008333333, 0.0008333333 (x, y), extent : 90, 95, 20, 25 (xmin, xmax, ymin, ymax), crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

source : C:/Users/hafez/personal/GIS DATA/srtm_55_08/srtm_55_08/srtm_55_08.tif, names : srtm_55_08, values : -32768, 32767 (min, max)

hillshade

hill<-hillShade(slope,aspect,40,270

plot(hill,col = grey(0:100/100), legend = FALSE, main = “Hillshade”)
plot(r1, col = rainbow(25, alpha = 0.35), add = TRUE)

2: I want to get latitude, longitude, and elevation value from this file. code is here, now you can save it in your local machine as you desired file.

# lat long
ed = as.data.frame(r1,xy=TRUE)
# remove na value
ed2 = ed[!is.na(ed[,3]),]

you are thinking about missing something ? yeah? 3D? it is possible in R ,lets’ draw 3D DEM

3D DEM

library(plot3D)

plot3D(r1,zfac=0.5, drape=NULL, col=terrain.colors(5),
at=10, rev=FALSE,
useLegend=TRUE,main=”Digital Elevation Model (DEM)”)

3: I want to get elevation data and plot from the DEM file with python. it is easy, it seems to be tedious because the scrips contain some loops and nested loop. be patient, it gonna okay after few hours. you should be careful because of python use indentation. don’t worry, I have attached all codes in the GitHub repo and you can get it.

There are many ways to get elevation data from DEM. I am going to get through using transects. I drew 24 transects over the DEM file and then I wrote scripts. It will automatically generate a CSV file containing elevation data in the meter.

In this figure, the red line indicates Transect over DEM.

code: this is so elaborative code in python

import rasterio
import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

line=gpd.read_file(‘polyline.shp’)

start_coords=list(line.geometry[0].coords)[0]
end_coords=list(line.geometry[0].coords)[1]

n_points=22

lon=start_coords[0]
lat=start_coords[1]

for i in np.arange(1,n_points):
x_dist=end_coords[0]-start_coords[0]
y_dist=end_coords[0]-start_coords[0]
point=

for ind, row in line.iterrows():

XS_ID = row[‘Id’]

start_coords = list([row.geometry][0].coords)[0]
end_coords = list([row.geometry][0].coords)[1]

lon = [start_coords[0]]
lat = [start_coords[1]]

n_points = 22

for i in np.arange(1, n_points+1):
x_dist = end_coords[0] — start_coords[0]
y_dist = end_coords[1] — start_coords[1]
point = [(start_coords[0] + (x_dist/(n_points+1))*i), (start_coords[1] + (y_dist/(n_points+1))*i)]
lon.append(point[0])
lat.append(point[1])

lon.append(end_coords[0])
lat.append(end_coords[1])
df = pd.DataFrame({‘Latitude’: lat,
‘Longitude’: lon})
gdf_pcs= gpd.GeoDataFrame(df, geometry = gpd.points_from_xy(df.Longitude, df.Latitude))
gdf.crs = {‘init’: ‘epsg:4682’}
#gdf_pcs = gdf.to_crs(epsg = 4682)

gdf_pcs[‘h_distance’] = 0
for index, row in gdf_pcs.iterrows():
gdf_pcs[‘h_distance’].loc[index] = gdf_pcs.geometry[0].distance(gdf_pcs.geometry[index])

# Extracting the elevations from the DEM

gdf_pcs[‘Elevation’] = 0
dem = rasterio.open(r’srtm_55_08.tif’, mode = ‘r’)
for index, row in gdf_pcs.iterrows():
row, col = dem.index(row[‘Longitude’], row[‘Latitude’])
dem_data = dem.read(1)
gdf_pcs[‘Elevation’].loc[index] = dem_data[row, col]
# Extract h_distance (x) and Elevation (y) columns into a Pandas DataFrame
x_y_data = gdf_pcs[[‘h_distance’, ‘Elevation’]]
x_y_data.to_csv(r’extracted_sections’+ str(XS_ID) + ‘.csv’ )
# Creating plots for each cross sectional profile
plt.plot(gdf_pcs[‘h_distance’], gdf_pcs[‘Elevation’])
plt.xlabel(‘Distance (m)’)
plt.ylabel(‘Elevation (m)’)
plt.grid(True)
plt.title(XS_ID)
plt.savefig(r’extracted_sections’ + str(XS_ID) + ‘.png’ )
plt.show()

do you gonna see what I get after this?

CSV file and profile plot

you wanna plot in python, let's start coding

general plotting

histogram

code:

f=rasterio.open(‘srtm_55_08.tif’)

from rasterio.plot import show

show(f)

from rasterio.plot import show_hist

show_hist(f,bins=50, lw=0.0, stacked=False, alpha=0.3, histtype=’stepfilled’, title=”Histogram”)

Do you wanna contour plot ? because we studied it a lot….

contour!!

code:

from osgeo import gdal
gdal_data = gdal.Open(‘srtm_55_08.tif’)
gdal_band = gdal_data.GetRasterBand(1)
nodataval = gdal_band.GetNoDataValue()

# convert to a numpy array
data_array = gdal_data.ReadAsArray().astype(np.float)
data_array

# replace missing values if necessary
if np.any(data_array == nodataval):
data_array[data_array == nodataval] = np.nan
#Plot out data with Matplotlib’s ‘contour’
fig = plt.figure(figsize = (12, 8))
ax = fig.add_subplot(111)
plt.contour(data_array, cmap = “terrain”,
levels = list(range(0, 5000, 100))) # you can change this value , explore it
plt.title(“Elevation Contours of Chittagong “)
cbar = plt.colorbar()
plt.gca().set_aspect(‘equal’, adjustable=’box’)
plt.show()

#contourf’ to see filled contours.

fig = plt.figure(figsize = (12, 8))
ax = fig.add_subplot(111)
plt.contourf(data_array, cmap = “terrain”,
levels = list(range(0, 5000, 100)))
plt.title(“Elevation Contours of Chittagong [fill]”)
cbar = plt.colorbar()
plt.gca().set_aspect(‘equal’, adjustable=’box’)
plt.show()

Okay,guys, I wanna do the same thing what I did with R . I am talking about hillside mapping. okay, no more talk,

bounds=f.bounds

extent = xmin, xmax, ymin, ymax = 90.0, 95.0, 20.0, 25.0

def hillshade(array, azimuth, angle_altitude):

x, y = np.gradient(array)
slope = np.pi/2. — np.arctan(np.sqrt(x*x + y*y))
aspect = np.arctan2(-x, y)
azimuthrad = azimuth*np.pi / 180.
altituderad = angle_altitude*np.pi / 180.

shaded = np.sin(altituderad) * np.sin(slope) \
+ np.cos(altituderad) * np.cos(slope) \
* np.cos(azimuthrad — aspect)
return 255*(shaded + 1)/2
topocmap = plt.get_cmap(‘Spectral_r’)
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111)
ax.matshow(hillshade(data_array, 30, 30), extent=extent, cmap=’Greys’, alpha=.5, zorder=10)
cax = ax.contourf(data_array, np.arange(vmin, vmax, 10),extent=extent,
cmap=topocmap, vmin=vmin, vmax=vmax, origin=’image’)
fig.colorbar(cax, ax=ax)
fig.savefig(‘images.png’)

I have attached Python [version 3.7.0] code in the GitHub repository. Please check it. The next post will be [1. Practical Ecological data analysis

2. Cluster and principal component analysis

3. Analysis of biodiversity data

4. The cross-sectional plot of ocean

5. Temperature extraction [in CSV format], time series generation and visualization from NetCDF file]

6: species distribution modeling

If you have any question, please feel free to ask. And if you can improve the code for me, please do. I really appreciate it. Thank you.