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.

In [ ]:

```
# 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.

In [ ]:

```
# 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.

In [ ]:

```
ds_dem.plot()
```

In [ ]:

```
ds_dem.values
```

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

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

```
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

In [ ]:

```
# Select only data that are below our chosen point
data_point_down = data_point.where(ds_dem < elevation_point)
```

In [ ]:

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

In [ ]:

```
# Select only data that are above our chosen point
data_point_up = data_point.where(ds_dem > elevation_point)
```

In [ ]:

```
# 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.

First we load the raion shapefile available in the `data`

folder

In [ ]:

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

In [ ]:

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

In [ ]:

```
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`

In [ ]:

```
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.)

In [ ]:

```
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.

In [ ]:

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