Alvin C.G. Varquez (varquez.a.aa@m.titech.ac.jp)
Here, we will try to explore the ways to accomplish the following tasks.
Global_Landslide_Catalog_Export.csv
.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:
rasterio
and geopandas
.!
) 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.
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.
!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.
Here, we will look at common ways to visualize a raster and shapefile in Google Colab.
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.
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.
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).
!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)
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).
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.
geo_df.plot()
<Axes: >
The following is how to focus only on a certain region.
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
(-25.0, 50.0)
The following utilizes a module called contextily
to add a basemap underneath the plot.
!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)
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)
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.
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])
(-23.0, -22.0)
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:3857
italicized text
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.
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
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
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).
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.
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.
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.
!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.
!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.
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.
!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
To calculate the slope and hillshade from the projected DEM/DSM file, we can use the following commands separately.
!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.
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])
(-2632033.410913339, -2511525.2348457137)
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])
(-2632033.410913339, -2511525.2348457137)
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.
!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.
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.
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
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.
!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%)
!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
.
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.
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)
# "out" can be saved into another shapefile.
out[out['count']>0]