More Advanced Data Cube Usage

Introduction on how to manipulate DEMs

In this notebook you will learn how to manipulate Digital Elevation Models (DEMs) using python.

This notebook will require a number of external Python modules to run, which will be imported at the top of the notebook as his convention.

# Load all the packages needed
import sys
import geopandas
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
sys.path.append('/eOsphere/bin')
from dcFunctions import coord_convert

Let's load the DEM DEM_KG_UTM43.tif available in the data folder. Just for the sake of this notebook, this DEM has been reprojected to UTM 43N.

# Load the reprojected DEM with xarray

ds_dem = xr.open_rasterio('./data/DEM_KG_UTM43.tif')
print(ds_dem)

If we try to plot the data at this stage we will find that we have to mask the nodata value, to see the elevation of the terrain.

ds_dem.plot()

ds_dem.values

For analysis and plotting we are going to remove the nodata value, so that it will appear as NaN

# Find the nodata value to remove it from our DEM data
nodata = ds_dem.nodatavals[0]
print(nodata)
print()

# Remove the nodata value
ds_dem = ds_dem.where(ds_dem!=nodata)
print(ds_dem)

# Plotting DEM
ds_dem.plot(figsize=(18,8))

Now you can treat the DEM like a data layer, so for example you could find the maximum and minimum elevation point for Kyrgyzstan

print('Minimum elevation point:')
print(ds_dem.min().values, 'm')
print()
print('Maximum elevation point:')
print(ds_dem.max().values, 'm')

Find a region around a point

# Convert a chosen point in UTM43 coordinates
lon, lat = 77, 41.7
x_pt, y_pt = coord_convert(lon, lat, 4326, 32643)
x_pt, y_pt

# Plot of the reprojected DEM with the chosen point 
# marked with an orange cross 

ds_dem.plot(figsize=(18,8))
plt.plot(x_pt,y_pt,'+r',ms=20)

Now we would like to select an area around the selected point, for example an area around the point with a radius of 50 km

# Select only the data within a 50 km radium from our chosen point

r_coverage = 50000. # 50 km radius
# Circle: (x-x0)^2 + (y-y0)^2 = r^2
condition_area = (ds_dem.x - x_pt)**2. + (ds_dem.y - y_pt)**2 < r_coverage**2

data_point = ds_dem.where(condition_area, drop=True)

# Plot of our selected data
data_point.plot(figsize=(11,9))
plt.plot(x_pt,y_pt,'+r',ms=10)

Let's find the elevation of the chosen point

elevation_point = (ds_dem.sel(x=x_pt, method='nearest')).sel(y=y_pt, method='nearest') 
print(elevation_point.values)

Now, we can plot only the points that are below or about the chosen point

# Select only data that are below our chosen point 

data_point_down = data_point.where(ds_dem < elevation_point)

# Plot data below our point within 50 km 

data_point_down.plot(figsize=(11,9), vmin = 1200, vmax=4700)
plt.plot(x_pt,y_pt,'+r',ms=10)

# Select only data that are above our chosen point 

data_point_up = data_point.where(ds_dem > elevation_point)

# Plot data above our point within 50 km 

data_point_up.plot(figsize=(11,9), vmin = 1200, vmax=4700)
plt.plot(x_pt,y_pt,'+r',ms=10)

This selection of data can have several applications. For example, if we have the position of the meteo stations in Kyrgyzstan, we could check which part of the country is covered.

Now, we are going to see how to select the data if we have multiple points.

Given a list of points, let's find a region around each single point

First we load the raion shapefile available in the data folder

# Loading again the Raion shapefiles
raion_gdf = geopandas.read_file("data/kgz_admbnda_adm2_moes_20181119.shp", encoding='utf-8')
raion_utm_gdf = raion_gdf.to_crs({'init': 'epsg:32643'})

Then, we create a list of points and convert the coordinate in UTM43N (you can also load a shapefile of points).

longs = [70.8, 71.5, 73.5, 74, 75, 76, 76.7, 78]
lats = [40, 42.5, 42, 41, 41.6, 42.5, 41.2, 42]
x_pts, y_pts = np.zeros(len(longs)), np.zeros(len(longs))
for i in range(len(longs)):
    x_pts[i], y_pts[i] = coord_convert(longs[i], lats[i], 4326, 32643)

ds_dem.plot(figsize=(18,8))
for i in range(len(longs)):
    plt.plot(x_pts[i],y_pts[i], color='red', marker='o', markersize=20, markerfacecolor='none')

Now, we want to apply the condition used in the previous cell, where we select only a circle around the point, for all the points.

We can start with a first point to initialise our condition_all_pt

x_pt, y_pt = x_pts[0], y_pts[0]
condition_all_pt = (ds_dem.x - x_pt)**2. + (ds_dem.y - y_pt)**2 < r_coverage**2

Once we have our first point, we can loop for all the points (x_pts and y_pts) and sum all the conditions together. (Summing them we can keep all the points.)

for i in range(1,len(x_pts)):
    x_pt, y_pt = x_pts[i], y_pts[i]
    condition_temp = (ds_dem.x - x_pt)**2. + (ds_dem.y - y_pt)**2 < r_coverage**2
    condition_all_pt += condition_temp
data_stations = ds_dem.where(condition_all_pt)

From the plot we can see that which parts of Kygyzstan are not being covered.

fig, ax = plt.subplots(figsize=(18,8))
data_stations.plot(ax=ax)
raion_utm_gdf.plot(ax=ax, facecolor='none', edgecolor='black')
for i in range(len(longs)):
    plt.plot(x_pts[i],y_pts[i], color='red', marker='+', markersize=10)