Lecture 4 Supplementary Handouts¶

Alvin C.G. Varquez (varquez.a.aa@m.titech.ac.jp)

I. Introduction¶

Here, we will try to explore the ways to accomplish the following tasks.

  1. Select a region that experienced more than 10 landslides. Download the "comma-separated values file" here. We will use the downloaded csv file Global_Landslide_Catalog_Export.csv.
  2. Download the DEM tile from JAXA that covers the region. See in instructions in "4. Download" inside the link.
  3. Create a 50-m radius buffer around the landslide events, calculate and summarize the statistics of the slopes contained in all buffers.
  4. Based on the statistics and visualization of data, discuss the possible risky slope for landslides to occur and the look at other factors that can trigger landslides. Include news clippings of past events.
  5. Propose a strategy to mitigate the incidents of landslides by citing existing strategies.

We will use Google Colab notebooks to conduct this analyses. Note that you can also install this system into your computer but I chose to not request this to avoid OS compatibility issues.

In order to meet the tasks above, it is necessary to learn the following:

  • Basic visualization of the raw data with rasterio and geopandas.
  • Creating vector from comma-separated values (this can also be done in QGIS)
  • Converting the projection
  • Using the GDAL toolbox programatically (in both Python and in the background using !) that is the underlying toolbox of QGIS.

The whole exercise can be summed up by a rough-sketch I made. Make sure to upload all the necessary files into the Google Colab session. Also, as you'll notice in each of the cells, they are designed to run on its own. In some cells, the modules are also called. This means that the whole cell can be saved as a separate python script, given, of course, that the modules were installed in advance. pip may be used to install these modules.

Flowchart

Installing all necessary modules¶

Each time you open your notebook on Google Colab, your system connects to a remote server that does not necessarily have the modules installed by default. Hence, it is necessary to install them in the background. For example, the following installs a module named rasterio in the remote server. A command that starts with ! is a unix command that is run on the remote server and not in the python script.

In [1]:
!pip install rasterio ## Installs the module rasterio
!apt install gdal-bin python3-gdal --quiet ## Installs the module for gdal

Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Requirement already satisfied: rasterio in /usr/local/lib/python3.10/dist-packages (1.3.6)
Requirement already satisfied: click-plugins in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.1.1)
Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from rasterio) (67.7.2)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from rasterio) (2022.12.7)
Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.10/dist-packages (from rasterio) (8.1.3)
Requirement already satisfied: affine in /usr/local/lib/python3.10/dist-packages (from rasterio) (2.4.0)
Requirement already satisfied: snuggs>=1.4.1 in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.4.7)
Requirement already satisfied: numpy>=1.18 in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.22.4)
Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.10/dist-packages (from rasterio) (0.7.2)
Requirement already satisfied: attrs in /usr/local/lib/python3.10/dist-packages (from rasterio) (23.1.0)
Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.10/dist-packages (from snuggs>=1.4.1->rasterio) (3.0.9)
Reading package lists...
Building dependency tree...
Reading state information...
gdal-bin is already the newest version (3.3.2+dfsg-2~focal2).
python3-gdal is already the newest version (3.3.2+dfsg-2~focal2).
0 upgraded, 0 newly installed, 0 to remove and 26 not upgraded.

II. Visualizing and Creating Shapefiles¶

Here, we will look at common ways to visualize a raster and shapefile in Google Colab.

Visualizing Raster¶

To visualize a raster, the following code may be used that uses the rasterio module that was installed earlier. I have uploaded the tif file into Google Colab and copied the path (just drag the file into the folder on the left). Note that it is also possible to link the notebook to your Google drive folder.

In [2]:
import rasterio
import rasterio.plot

fp = r'/content/ALPSMLC30_S023W044_DSM.tif'
img = rasterio.open(fp)
ax = rasterio.plot.show(img, vmin=0.,vmax=800.0,cmap='jet')

The image above shows the tile's DSM with the axes representing longitude and latitude. Tip: In order to print the boundaries of the tif above, insert print(img.bounds) in the code.

Visualizing Vectors and creating a shapefile from the csv file.¶

To visualize vectors and create vectors from text files, such as the Global_Landslide_Catalog_Export.csv.csv file, we can use a module called geopandas. Again, they can be installed for use into Google Colab with the following command (use this only once).

In [3]:
!pip install geopandas

Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Requirement already satisfied: geopandas in /usr/local/lib/python3.10/dist-packages (0.13.0)
Requirement already satisfied: pandas>=1.1.0 in /usr/local/lib/python3.10/dist-packages (from geopandas) (1.5.3)
Requirement already satisfied: packaging in /usr/local/lib/python3.10/dist-packages (from geopandas) (23.1)
Requirement already satisfied: shapely>=1.7.1 in /usr/local/lib/python3.10/dist-packages (from geopandas) (2.0.1)
Requirement already satisfied: fiona>=1.8.19 in /usr/local/lib/python3.10/dist-packages (from geopandas) (1.8.22)
Requirement already satisfied: pyproj>=3.0.1 in /usr/local/lib/python3.10/dist-packages (from geopandas) (3.5.0)
Requirement already satisfied: click-plugins>=1.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (1.1.1)
Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (67.7.2)
Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (8.1.3)
Requirement already satisfied: attrs>=17 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (23.1.0)
Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (0.7.2)
Requirement already satisfied: munch in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (2.5.0)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (2022.12.7)
Requirement already satisfied: six>=1.7 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (1.16.0)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (2022.7.1)
Requirement already satisfied: numpy>=1.21.0 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (1.22.4)
Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (2.8.2)

Converting the text file into a shapefile.¶

In the program below, we used pandas module (aliased as pd) to load the csv file as a dataframe object df. This dataframe contains a column of latitude and longitude. They can be confirmed separately by printing the columns of the dataframe using print(df.columns).

One thing that makes vector file unique is that it contains a geometry information. It is necessary to derive a geometry from the latitude and longitude pairs. This is achieved with the Point module and the line where geometry = ... is located. Here, the latitude and longitude pairs are used to obtain geometry. Once the geometry is estimated, a new dataframe geo_df that contains both geometry and the projection (for this case EPSG:4326) is created. geo_df is essentially a pandas-type of dataframe but with a projection and geometry. The selection of the projection was based on intuition that the latitude and longitude units are in arc-degrees.

geo_df.plot() is used to draw the points over the map. geo_df.to_file('output.shp') may be used to save the current dataframe into a shapefile named landslide_wgs84.shp (can be visualized or processed using QGIS as well).

In [4]:
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point

df = pd.read_csv('/content/Global_Landslide_Catalog_Export.csv')

#The following loops through each latitude and longitude to create a geometry column.
geometry = [Point(xy) for xy in zip(df['longitude'],df['latitude'])]

#This creates a new dataframe with the geometry column calculated above.
geo_df = gpd.GeoDataFrame(df,crs='epsg:4326',geometry=geometry)

geo_df.to_file('landslide_wgs84.shp', encoding='utf-8')

<ipython-input-4-53c3f3880b2a>:13: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.
  geo_df.to_file('landslide_wgs84.shp', encoding='utf-8')
WARNING:fiona._env:Normalized/laundered field name: 'source_name' to 'source_nam'
WARNING:fiona._env:Normalized/laundered field name: 'source_link' to 'source_lin'
WARNING:fiona._env:Normalized/laundered field name: 'event_title' to 'event_titl'
WARNING:fiona._env:Normalized/laundered field name: 'event_description' to 'event_desc'
WARNING:fiona._env:Normalized/laundered field name: 'location_description' to 'location_d'
WARNING:fiona._env:Normalized/laundered field name: 'location_accuracy' to 'location_a'
WARNING:fiona._env:Normalized/laundered field name: 'landslide_category' to 'landslide_'
WARNING:fiona._env:Normalized/laundered field name: 'landslide_trigger' to 'landslid_1'
WARNING:fiona._env:Normalized/laundered field name: 'landslide_size' to 'landslid_2'
WARNING:fiona._env:Normalized/laundered field name: 'landslide_setting' to 'landslid_3'
WARNING:fiona._env:Normalized/laundered field name: 'fatality_count' to 'fatality_c'
WARNING:fiona._env:Normalized/laundered field name: 'injury_count' to 'injury_cou'
WARNING:fiona._env:Normalized/laundered field name: 'event_import_source' to 'event_impo'
WARNING:fiona._env:Normalized/laundered field name: 'event_import_id' to 'event_im_1'
WARNING:fiona._env:Normalized/laundered field name: 'country_name' to 'country_na'
WARNING:fiona._env:Normalized/laundered field name: 'country_code' to 'country_co'
WARNING:fiona._env:Normalized/laundered field name: 'admin_division_name' to 'admin_divi'
WARNING:fiona._env:Normalized/laundered field name: 'admin_division_population' to 'admin_di_1'
WARNING:fiona._env:Normalized/laundered field name: 'gazeteer_closest_point' to 'gazeteer_c'
WARNING:fiona._env:Normalized/laundered field name: 'gazeteer_distance' to 'gazeteer_d'
WARNING:fiona._env:Normalized/laundered field name: 'submitted_date' to 'submitted_'
WARNING:fiona._env:Normalized/laundered field name: 'created_date' to 'created_da'
WARNING:fiona._env:Normalized/laundered field name: 'last_edited_date' to 'last_edite'
WARNING:fiona._env:Value 'Another landslide in sitio Bakilid in barangay Lahug also left two families homeless. Lilibeth Magsuling was breastfeeding her two-month-old baby outside their house at 12 noon yesterday when she heard a loud blast behind her. She immediately hugged her baby boy and ran. When she looked back, she saw bamboo trees and soil covering her house. “It’s fortunate that we were outside the house. I was about to put my baby in the cradle inside the house,” she said. Elenit Villaflor was also outside her house in sitio Bakilid with her two-month-old baby when the landslide struck. ' of field event_desc has been truncated to 254 characters.  This warning will not be emitted any more for that layer.

Note that the warnings above mean that there are column names that were too long and were shortened.

To visualize the data, the following command will do.

In [5]:
geo_df.plot()
Out[5]:
<Axes: >

The following is how to focus only on a certain region.

In [6]:
import matplotlib.pyplot as plt
geo_df.plot()
plt.xlim([100.0,150.0]) #Range of longitude
plt.ylim([-25.0,50.0]) #Range of latitude
Out[6]:
(-25.0, 50.0)

The following utilizes a module called contextily to add a basemap underneath the plot.

In [7]:
!pip install contextily

Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Requirement already satisfied: contextily in /usr/local/lib/python3.10/dist-packages (1.3.0)
Requirement already satisfied: rasterio in /usr/local/lib/python3.10/dist-packages (from contextily) (1.3.6)
Requirement already satisfied: pillow in /usr/local/lib/python3.10/dist-packages (from contextily) (8.4.0)
Requirement already satisfied: mercantile in /usr/local/lib/python3.10/dist-packages (from contextily) (1.2.1)
Requirement already satisfied: matplotlib in /usr/local/lib/python3.10/dist-packages (from contextily) (3.7.1)
Requirement already satisfied: xyzservices in /usr/local/lib/python3.10/dist-packages (from contextily) (2023.2.0)
Requirement already satisfied: joblib in /usr/local/lib/python3.10/dist-packages (from contextily) (1.2.0)
Requirement already satisfied: geopy in /usr/local/lib/python3.10/dist-packages (from contextily) (2.3.0)
Requirement already satisfied: requests in /usr/local/lib/python3.10/dist-packages (from contextily) (2.27.1)
Requirement already satisfied: geographiclib<3,>=1.52 in /usr/local/lib/python3.10/dist-packages (from geopy->contextily) (2.0)
Requirement already satisfied: numpy>=1.20 in /usr/local/lib/python3.10/dist-packages (from matplotlib->contextily) (1.22.4)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib->contextily) (23.1)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib->contextily) (2.8.2)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib->contextily) (1.4.4)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib->contextily) (1.0.7)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib->contextily) (3.0.9)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib->contextily) (0.11.0)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib->contextily) (4.39.3)
Requirement already satisfied: click>=3.0 in /usr/local/lib/python3.10/dist-packages (from mercantile->contextily) (8.1.3)
Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from rasterio->contextily) (67.7.2)
Requirement already satisfied: click-plugins in /usr/local/lib/python3.10/dist-packages (from rasterio->contextily) (1.1.1)
Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.10/dist-packages (from rasterio->contextily) (0.7.2)
Requirement already satisfied: affine in /usr/local/lib/python3.10/dist-packages (from rasterio->contextily) (2.4.0)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from rasterio->contextily) (2022.12.7)
Requirement already satisfied: attrs in /usr/local/lib/python3.10/dist-packages (from rasterio->contextily) (23.1.0)
Requirement already satisfied: snuggs>=1.4.1 in /usr/local/lib/python3.10/dist-packages (from rasterio->contextily) (1.4.7)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /usr/local/lib/python3.10/dist-packages (from requests->contextily) (1.26.15)
Requirement already satisfied: charset-normalizer~=2.0.0 in /usr/local/lib/python3.10/dist-packages (from requests->contextily) (2.0.12)
Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.10/dist-packages (from requests->contextily) (3.4)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.7->matplotlib->contextily) (1.16.0)
In [8]:
import contextily as cx
ax = geo_df.plot()
plt.xlim([135,142.0])
plt.ylim([34,40])
cx.add_basemap(ax, crs=geo_df.crs)

Overlaying a raster and a shapefile¶

The above methods of visualizing a raster and a shapefile may be combined with the help of the matplotlib module to display both the raster and shapefile. The following is the command that merges the two files together.

In [9]:
import rasterio
from rasterio.plot import show
import matplotlib.pyplot as plt
from shapely.geometry import Point
import geopandas as gpd
import pandas as pd

df = pd.read_csv('/content/Global_Landslide_Catalog_Export.csv')

geometry = [Point(xy) for xy in zip(df['longitude'],df['latitude'])]
geo_df = gpd.GeoDataFrame(df,crs='epsg:4326',geometry=geometry)

fp = r'/content/ALPSMLC30_S023W044_DSM.tif'
img = rasterio.open(fp)
fig,ax = plt.subplots(figsize=(10,12))

## Plotting the raster with a colorbar
ras = rasterio.plot.show(img, vmin=0.,vmax=1000.0,cmap='nipy_spectral',ax=ax)
im = ras.get_images()[0]
fig.colorbar(im, ax=ax,orientation='horizontal')

## Plotting the shapefile
geo_df.plot(facecolor='red', edgecolor='black',ax=ax)

## Constraining to the region of the raster
plt.xlim([img.bounds.left,img.bounds.right])
plt.ylim([img.bounds.bottom,img.bounds.top])
Out[9]:
(-23.0, -22.0)

III. Changing the projection of the shapefile from geographic to projected coordinate system.¶

Many web applications such as Google, Bing, ArcGIS use projected coordinate systems because it provide a more useful/practical units of distance. Specifically, EPSG:3857 is used.

We can copy the geo_df earlier to a dataframe geo_proj where we'll convert the EPSG:4326 to EPSG:3857italicized text

In [10]:
geo_proj = geo_df.copy()
geo_proj['geometry'] = geo_df['geometry'].to_crs('epsg:3857')
geo_proj.crs = 'epsg:3857'

If you print the geo_proj['geometry'] and geo_df['geometry'], you will notice the differences in value corresponding to the differences in the units.

In [11]:
print(geo_proj['geometry'])
print(geo_df['geometry'])

0         POINT (11961279.286 3837376.100)
1        POINT (-13654782.699 5687886.023)
2        POINT (-8388892.111 -1246795.733)
3          POINT (9095692.954 3355017.763)
4         POINT (13792240.006 1156618.626)
                       ...                
11028     POINT (12432154.542 -876987.936)
11029      POINT (8424727.107 3948929.094)
11030     POINT (10216016.984 3021591.066)
11031      POINT (8178907.765 4995599.361)
11032      POINT (8722606.264 1975095.120)
Name: geometry, Length: 11033, dtype: geometry
0         POINT (107.45000 32.56250)
1        POINT (-122.66300 45.42000)
2        POINT (-75.35870 -11.12950)
3          POINT (81.70800 28.83780)
4         POINT (123.89780 10.33360)
                    ...             
11028     POINT (111.67994 -7.85341)
11029      POINT (75.68061 33.40308)
11030      POINT (91.77204 26.18161)
11031      POINT (73.47238 40.88639)
11032      POINT (78.35651 17.46563)
Name: geometry, Length: 11033, dtype: geometry

IV. Creating a buffer geometry centered at the points from the point shapefile or geopandas dataframe.¶

To create a buffer geometry surrounding the points, we can revise the geometry of geo_proj by using the buffer function. Note that the value inside the buffer corresponds to the units of the projection of the geometry. This means that to create 50-m. buffers around the points, we can use the following command. (warning: be careful when running the cell below. Repeatedly using it might continually make buffers)link text

In [12]:
geo_proj['geometry'] = geo_proj.buffer(50.0)

We can save the shapefile containing buffers using the following command. If we use the shapefile below in the visualization with the raster above, the two files won't overlap in the same map because the DSM above uses a different coordinate reference system (CRS).

In [13]:
geo_proj.to_file('landslide_epsg3857.shp', encoding='utf-8')

<ipython-input-13-0f83e96d7145>:1: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.
  geo_proj.to_file('landslide_epsg3857.shp', encoding='utf-8')
WARNING:fiona._env:Normalized/laundered field name: 'source_name' to 'source_nam'
WARNING:fiona._env:Normalized/laundered field name: 'source_link' to 'source_lin'
WARNING:fiona._env:Normalized/laundered field name: 'event_title' to 'event_titl'
WARNING:fiona._env:Normalized/laundered field name: 'event_description' to 'event_desc'
WARNING:fiona._env:Normalized/laundered field name: 'location_description' to 'location_d'
WARNING:fiona._env:Normalized/laundered field name: 'location_accuracy' to 'location_a'
WARNING:fiona._env:Normalized/laundered field name: 'landslide_category' to 'landslide_'
WARNING:fiona._env:Normalized/laundered field name: 'landslide_trigger' to 'landslid_1'
WARNING:fiona._env:Normalized/laundered field name: 'landslide_size' to 'landslid_2'
WARNING:fiona._env:Normalized/laundered field name: 'landslide_setting' to 'landslid_3'
WARNING:fiona._env:Normalized/laundered field name: 'fatality_count' to 'fatality_c'
WARNING:fiona._env:Normalized/laundered field name: 'injury_count' to 'injury_cou'
WARNING:fiona._env:Normalized/laundered field name: 'event_import_source' to 'event_impo'
WARNING:fiona._env:Normalized/laundered field name: 'event_import_id' to 'event_im_1'
WARNING:fiona._env:Normalized/laundered field name: 'country_name' to 'country_na'
WARNING:fiona._env:Normalized/laundered field name: 'country_code' to 'country_co'
WARNING:fiona._env:Normalized/laundered field name: 'admin_division_name' to 'admin_divi'
WARNING:fiona._env:Normalized/laundered field name: 'admin_division_population' to 'admin_di_1'
WARNING:fiona._env:Normalized/laundered field name: 'gazeteer_closest_point' to 'gazeteer_c'
WARNING:fiona._env:Normalized/laundered field name: 'gazeteer_distance' to 'gazeteer_d'
WARNING:fiona._env:Normalized/laundered field name: 'submitted_date' to 'submitted_'
WARNING:fiona._env:Normalized/laundered field name: 'created_date' to 'created_da'
WARNING:fiona._env:Normalized/laundered field name: 'last_edited_date' to 'last_edite'
WARNING:fiona._env:Value 'Another landslide in sitio Bakilid in barangay Lahug also left two families homeless. Lilibeth Magsuling was breastfeeding her two-month-old baby outside their house at 12 noon yesterday when she heard a loud blast behind her. She immediately hugged her baby boy and ran. When she looked back, she saw bamboo trees and soil covering her house. “It’s fortunate that we were outside the house. I was about to put my baby in the cradle inside the house,” she said. Elenit Villaflor was also outside her house in sitio Bakilid with her two-month-old baby when the landslide struck. ' of field event_desc has been truncated to 254 characters.  This warning will not be emitted any more for that layer.

V. Processing geotiffs with the GDAL toolbox.¶

In this exercise, we will use the built-in GDAL commands to process the DSM geotiff. Specifically, we will calculate the hillshade and the slope.

Reprojecting the geotiff¶

When we look at locations of past landslides, we find that it's not actually the value of topography itself but rather the slope or gradient of the topographic data. To calculate the slope, it is better to use the metric units. Hence, we should first project the DSM geotiff. To do so, we can use the gdalwarp command. Note that this is not inside a Python script (since the line is preceeded by the !). This means we will run this command directly inside the server or cloud where we're running the notebook.

In [14]:
!gdalwarp

Usage: gdalwarp [--help-general] [--formats]
    [-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"]* [-novshiftgrid]
    [-order n | -tps | -rpc | -geoloc] [-et err_threshold]
    [-refine_gcps tolerance [minimum_gcps]]
    [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]
    [-ovr level|AUTO|AUTO-n|NONE] [-wo "NAME=VALUE"] [-ot Byte/Int16/...] [-wt Byte/Int16]
    [-srcnodata "value [value...]"] [-dstnodata "value [value...]"] -dstalpha
    [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]
    [-cutline datasource] [-cl layer] [-cwhere expression]
    [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]
    [-if format]* [-of format] [-co "NAME=VALUE"]* [-overwrite]
    [-nomd] [-cvmd meta_conflict_value] [-setci] [-oo NAME=VALUE]*
    [-doo NAME=VALUE]*
    srcfile* dstfile

Available resampling methods:
    near (default), bilinear, cubic, cubicspline, lanczos, average, rms,
    mode,  max, min, med, Q1, Q3, sum.

FAILURE: No target filename specified.

To convert the projection from EPGS:4326 of the raster file to EPSG:3857 we can use the following command.

In [15]:
!gdalwarp -s_srs 'epsg:4326' -t_srs 'epsg:3857' /content/ALPSMLC30_S023W044_DSM.tif projected_dem.tif

Processing /content/ALPSMLC30_S023W044_DSM.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

After running the program above, a raster file called projected_dem.tif is created that is essentially the DSM file but of a projected coordinate system.

Showing the details of the geotiff.¶

To show important details of the geotiff, we can use the following command. The output will show very useful information such as the size of the raster, the projection, the pixel size, and the boundaries.

In [16]:
!gdalinfo /content/projected_dem.tif

Driver: GTiff/GeoTIFF
Files: /content/projected_dem.tif
Size is 3455, 3740
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06°S and 85.06°N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-4898057.594904037192464,-2511525.234845713712275)
Pixel Size = (32.221437451236639,-32.221437451236639)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-4898057.595,-2511525.235) ( 44d 0' 0.00"W, 22d 0' 0.00"S)
Lower Left  (-4898057.595,-2632033.411) ( 44d 0' 0.00"W, 23d 0' 0.44"S)
Upper Right (-4786732.529,-2511525.235) ( 42d59'59.82"W, 22d 0' 0.00"S)
Lower Right (-4786732.529,-2632033.411) ( 42d59'59.82"W, 23d 0' 0.44"S)
Center      (-4842395.062,-2571779.323) ( 43d29'59.91"W, 22d30' 3.47"S)
Band 1 Block=3455x1 Type=Int16, ColorInterp=Gray

Calculating the slope and hillshade¶

To calculate the slope and hillshade from the projected DEM/DSM file, we can use the following commands separately.

In [17]:
!gdaldem slope -p projected_dem.tif slope.tif
!gdaldem hillshade projected_dem.tif hillshade.tif

0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.

The commands above will calculate the slope and hillshade of the projected DEM/DSM file. They can be visualized by the following cells.

In [18]:
import rasterio
from rasterio.plot import show
import matplotlib.pyplot as plt
import geopandas as gpd

geo_df = gpd.read_file('/content/landslide_epsg3857.shp')

fp = r'/content/hillshade.tif'
img = rasterio.open(fp)
fig,ax = plt.subplots(figsize=(10,10))
rasterio.plot.show(img, vmin=0.,vmax=800.0,cmap='Reds',ax=ax)
geo_df.plot(facecolor='red', edgecolor='black',ax=ax)
plt.xlim([img.bounds.left,img.bounds.right])
plt.ylim([img.bounds.bottom,img.bounds.top])
Out[18]:
(-2632033.410913339, -2511525.2348457137)
In [19]:
import rasterio
from rasterio.plot import show
import matplotlib.pyplot as plt
import geopandas as gpd

geo_df = gpd.read_file('/content/landslide_epsg3857.shp')

fp = r'/content/slope.tif'
img = rasterio.open(fp)
fig,ax = plt.subplots(figsize=(10,10))
rasterio.plot.show(img, vmin=0.,vmax=100.0,cmap='nipy_spectral',ax=ax)
geo_df.plot(facecolor='red', edgecolor='red',ax=ax)
plt.xlim([img.bounds.left,img.bounds.right])
plt.ylim([img.bounds.bottom,img.bounds.top])
Out[19]:
(-2632033.410913339, -2511525.2348457137)

Extracting slope values that are contained inside the buffer¶

A convenient method to achieve this is again found in GDAL library called gdal_rasterize. This command can edit an existing geotiff using a polygon. The following command will replace the values of pixels not contained in the shapefile with a value of -99. Note that we first copy a slope file into another file (this copied file will be used). If we want the opposite where we replace the values of pixels bounded by the buffers, we simply remove -i in the command.

In [20]:
!cp /content/slope.tif ./slope_burn.tif
!gdal_rasterize -i -burn -99 -l landslide_epsg3857 landslide_epsg3857.shp /content/slope_burn.tif

0...10...20...30...40...50...60...70...80...90...100 - done.

The slope_burn.tif will be further analyzed to calculate the statistics of the slopes where the landslides in the past have occurred.

Reading the values of the geotiffs into an array for calculating the statistics.¶

We can use rasterio to read the contents of the slope_burn.tif into an array. We can then use common python libraries such as numpy to estimate the statistics of the pixel values. In the process, we also exclude the pixels that contain the value of -99 in the statistics (since pixels with -99 values are not within the buffers).

The opened raster file contained in the src contains the value information of each pixel. They can be stored into an array using the read function with 1 meaning to read the first band of the raster.

In [21]:
import rasterio
import numpy as np

src = rasterio.open('/content/slope_burn.tif')
array = src.read(1)
samples = array[array!=-99]
print('Mean:',np.mean(samples))
print('75th percentile:',np.percentile(samples,75))
print('25th percentile:',np.percentile(samples,25))
print('Median:',np.percentile(samples,50))

Mean: -1.9947972
75th percentile: 42.809889793395996
25th percentile: 13.743168115615845
Median: 29.05670928955078

VI. BONUS: Compressing all files for download¶

Sometimes downloading the files in Google colab one-by-one is not convenient. One way to do this is to compress all the files contained in the /content/ folder.

This can be achieved by running the zip command. The files in the /content/ folder will be stored in the outputfiles.zip, which can be downloaded and then processed locally in your computer.

In [22]:
!zip outputfiles.zip /content/*

updating: content/ALPSMLC30_S023W044_DSM.tif (deflated 48%)
updating: content/Global_Landslide_Catalog_Export.csv (deflated 64%)
updating: content/hillshade.tif (deflated 16%)
updating: content/landslide_epsg3857.cpg (stored 0%)
updating: content/landslide_epsg3857.dbf (deflated 91%)
updating: content/landslide_epsg3857.prj (deflated 42%)
updating: content/landslide_epsg3857.shp (deflated 54%)
updating: content/landslide_epsg3857.shx (deflated 63%)
updating: content/landslide_wgs84.cpg (stored 0%)
updating: content/landslide_wgs84.dbf (deflated 91%)
updating: content/landslide_wgs84.prj (deflated 17%)
updating: content/landslide_wgs84.shp (deflated 52%)
updating: content/landslide_wgs84.shx (deflated 71%)
updating: content/projected_dem.tif (deflated 49%)
updating: content/sample_data/ (stored 0%)
updating: content/slope_burn.tif (deflated 100%)
updating: content/slope.tif (deflated 49%)
  adding: content/drive/ (stored 0%)

BONUS: Conducting Zonal statistics in Python and not in QGIS¶

In [23]:
!pip install rasterstats

Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Requirement already satisfied: rasterstats in /usr/local/lib/python3.10/dist-packages (0.18.0)
Requirement already satisfied: numpy>=1.9 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (1.22.4)
Requirement already satisfied: cligj>=0.4 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (0.7.2)
Requirement already satisfied: simplejson in /usr/local/lib/python3.10/dist-packages (from rasterstats) (3.19.1)
Requirement already satisfied: fiona<1.9 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (1.8.22)
Requirement already satisfied: click>7.1 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (8.1.3)
Requirement already satisfied: rasterio>=1.0 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (1.3.6)
Requirement already satisfied: shapely in /usr/local/lib/python3.10/dist-packages (from rasterstats) (2.0.1)
Requirement already satisfied: affine<3.0 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (2.4.0)
Requirement already satisfied: attrs>=17 in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (23.1.0)
Requirement already satisfied: six>=1.7 in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (1.16.0)
Requirement already satisfied: click-plugins>=1.0 in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (1.1.1)
Requirement already satisfied: munch in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (2.5.0)
Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (67.7.2)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (2022.12.7)
Requirement already satisfied: snuggs>=1.4.1 in /usr/local/lib/python3.10/dist-packages (from rasterio>=1.0->rasterstats) (1.4.7)
Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.10/dist-packages (from snuggs>=1.4.1->rasterio>=1.0->rasterstats) (3.0.9)

The following will calculate the Zonal statistics and store it into a series called zonalstatistics_out.

In [ ]:
from rasterstats import zonal_stats
zonalstatistics_out = zonal_stats('/content/landslide_epsg3857.shp','/content/slope.tif')

We can then add the calculated zonalstatistics_out to each landslide point features by using the pd.concat function that merges according to the number of columns or rows.

In [ ]:
import geopandas as gpd
geo_landslide = gpd.read_file('/content/landslide_epsg3857.shp')
df = pd.DataFrame(zonalstatistics_out)
out = pd.concat([geo_landslide,df],axis=1)
In [ ]:
# "out" can be saved into another shapefile.
out[out['count']>0]