Skip to content

Commit

Permalink
Update intro_geoutils.py
Browse files Browse the repository at this point in the history
  • Loading branch information
@adehecq
adehecq committed Dec 13, 2023
1 parent 27235a0 commit 3290d1d
Showing 1 changed file with 166 additions and 48 deletions.
214 changes: 166 additions & 48 deletions notebooks/intro_geoutils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,72 +6,114 @@
# extension: .py
# format_name: light
# format_version: '1.5'
# jupytext_version: 1.14.0
# jupytext_version: 1.16.0
# kernelspec:
# display_name: Python 3 (ipykernel)
# language: python
# name: python3
# ---

# # Manipulating raster/vector data with geoutils
# # Manipulating raster/vector data in Python

# xdem relies on a second library, developed by the same group of people, to easily handle raster and vector data: **geoutils !**
# In this tutorial, we will learn some basic operations on raster/vector data in Python such as:
# - read, write, reproject and display rasters and vectors
# - calculate terrain attributes from DEMs (slope, aspect etc.)
# - calculate difference of DEMs and basic statistics.
#
# We will make use of the **geoutils** and **xdem** libraries. They are not "standard" Python libraries but were developed by glaciologists (R. Hugonnet, E. Mannerfelt, A. Dehecq), to make it easier to work with such data in Python. They are based upon the more known **rasterio/GDAL** and **geopandas** libraries.
# - **geoutils** provides functionalities for working with raster/vector data.
# - **xDEM** uses all functionnalities of geoutils and add additional tools specific to working with DEMs.
#
# The full documentation can be found at: \
# https://geoutils.readthedocs.io. \
# https://xdem.readthedocs.io.
#
# For more features, check the examples gallery on both geoutils and xdem's documentation.

# ### Import the necessary modules
# ## Import the necessary modules

# +
import matplotlib.pyplot as plt
import numpy as np
# %matplotlib widget

import geoutils as gu
import xdem
# -

# ### Download the sample data set (if not done already) - Should take a few seconds ###
# ## Download the sample data set (if not done already) - Should take a few seconds ###
# The code comes with some data for showcasing the functionality, here two DEMs over Longyearbyen, Norway and glacier outlines over Svalbard. \
# The files are already on disk, we only nned to find their location.

xdem.examples.download_longyearbyen_examples(overwrite=False)
print(xdem.examples.FILEPATHS_DATA["longyearbyen_ref_dem"])
print(xdem.examples.FILEPATHS_DATA["longyearbyen_tba_dem"])
print(xdem.examples.FILEPATHS_DATA["longyearbyen_glacier_outlines"])
print(xdem.examples.get_path("longyearbyen_ref_dem"))
print(xdem.examples.get_path("longyearbyen_tba_dem"))
print(xdem.examples.get_path("longyearbyen_glacier_outlines"))

# ### <span style='color:red '> **TO DO:** </span> Find these files on your computer and open them in QGIS.
#
#

# ## Read the two DEMs and glacier outlines for the region ###

# ### Read the two DEMs and glacier outlines for the region ###
# Raster files (e.g. GeoTiff) can be loaded in one line with `gu.Raster(path_to_file)`.

# Raster files (e.g. GeoTiff) can be loaded in one line with `gu.Raster(path_to_file)`. In xdem, the DEM class inherit the same functionalities and more, so here we use `xdem.DEM`.
dem_2009 = gu.Raster(xdem.examples.get_path("longyearbyen_ref_dem"))
dem_1990 = gu.Raster(xdem.examples.get_path("longyearbyen_tba_dem"))

dem_2009 = xdem.DEM(xdem.examples.FILEPATHS_DATA["longyearbyen_ref_dem"])
dem_1990 = xdem.DEM(xdem.examples.FILEPATHS_DATA["longyearbyen_tba_dem"])
# ### <span style='color:red '> **TO DO:** </span> Run the same commands, but copy/paste the file path. Uncomment the lines below (remove the # symbol) and replace ... with the file path.

# +
# dem_2009 = gu.Raster(...)
# dem_1990 = gu.Raster(...)
# -

# Vector files (e.g. ESRI shapefiles) can be loaded in one line with `gu.Vector(path_to_file)`.

outlines_1990 = gu.Vector(xdem.examples.FILEPATHS_DATA["longyearbyen_glacier_outlines"])
outlines_1990 = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))

# ### Quickly visualize a raster
# Since a Raster object comes with all atributes, it can be quickly plot with its georeferencing information.
# ## Quickly visualize a raster
# Since a Raster object comes with all atributes, it can be quickly plotted with its georeferencing information.

fig, ax = plt.subplots(figsize=(8, 6))
dem_2009.show()
plt.show()

# It is easier to visualize as a hillshade

# +
dem_2009_hs = xdem.terrain.hillshade(dem_2009)

fig, ax = plt.subplots(figsize=(8, 6))
dem_2009_hs.show(cmap='gray')
plt.show()
# -

# ### Quickly visualize vector data
# ## Quickly visualize vector data

outlines_1990.ds.plot()
fig, ax = plt.subplots(figsize=(8, 10))
outlines_1990.show(ax=ax)
plt.tight_layout()
plt.show()

# ### Notes on the Raster and Vector classes
# ## Notes on the Raster and Vector classes
#
# #### The `xdem.DEM` instances inherit from the `geoutils.Raster` class, which is based upon **rasterio**.
# ### The `geoutils.Raster` class is based upon **[rasterio](https://rasterio.readthedocs.io)** and inherits most of its metadata.
# These objects contain the raster metadata, with the same convention as rasterio.
# To georeference an object, one needs to know the coordinate reference system (called `crs` in rasterio/geopandas) and any 3 of these four information:
# - the raster width and height
# - the pixel resolution
# - the extent of the raster (called `bounds` in rasterio/geopandas)
# - the position of one pixel, traditionally, the upper-left
# - the position of one pixel, traditionally, the upper-left \
# These variables are inter-dependent, e.g. if one knows the raster's extent and width and height, the pixel resolution is fixed.
# All these variables are stored in the `xdem.DEM` instance with the following attributes:

print(dem_2009.crs.to_proj4()); print(dem_2009.crs.to_wkt())
# Coordinate reference system (CRS)

print(f"EPSG code: {dem_2009.crs.to_epsg()}\n")
print(f"Printed as WKT string: \n{dem_2009.crs.to_wkt()}")

# Raster width, height, resolution and bounds

dem_2009.width

Expand All @@ -87,82 +129,158 @@

# These information, and more, can all be obtained **at once** with the command

print(dem_2009)
print(dem_2009.info())

# or similarly with `dem_2009.info()`.
#
# Along with these metadata, the `xdem.DEM` object contains the data, stored as a numpy masked array in the `self.data` attribute:

dem_2009.data

# #### `gu.Vector` instances are based upon geopandas.
# The class contains several useful methods (`self.create_mask` is showcased below), and the GeoDataFrame can accessed via:
# ### `gu.Vector` instances are based upon geopandas.
# The class contains several useful methods (`self.create_mask` is showcased below), and the underlying GeoPandas' Data Frame can accessed via:

outlines_1990.ds

# ## Rasters operations
# # Terrain attributes
#
# **xDEM** enables calculating many terrain attributes: slope, aspect, hillshade, curbature etc. The full list is available here: https://xdem.readthedocs.io/en/stable/terrain.html
#
# Here we will demonstrate a few of them.

# ### Slope

slope = xdem.terrain.slope(dem_1990)
fig, ax = plt.subplots(figsize=(8, 6))
slope.show(cbar_title="Slope (degrees)")
plt.show()

# ### Aspect

aspect = xdem.terrain.aspect(dem_1990)
fig, ax = plt.subplots(figsize=(8, 6))
aspect.show(cbar_title="Aspect (degrees)", cmap="twilight")
plt.show()

# ### Reproject the two DEMs on the same grid
# #### <span style='color:red '> **TO DO:** </span>
# 1) Plot the terrain rugosity
# 2) Calculate the maximum slope (using `np.max`)
# Uncomment (remove the #) below.

# +
# rugosity = xdem.terrain.???(dem_1990)
# fig, ax = plt.subplots(figsize=(8, 6))
# rugosity.show(cbar_title="...", cmap="...")
# plt.show()

# +
# print("Maximum slope is:", np.max(...))
# -

# # Raster operations

# ## Reproject the two DEMs on the same grid
# NB: Here the function returns a warning because the two DEMs are already on the same grid, so nothing is actually done.

dem_1990 = dem_1990.reproject(dst_ref=dem_2009)
print(dem_1990)
print(dem_2009)

# ### Reproject to a given resolution, bounds, or CRS
# Check that both DEMs indeed have the same grid

print(dem_1990.info())
print(dem_2009.info())

# ## Reproject to a given resolution, bounds, or CRS

# #### Change pixel resolution

dem_test = dem_1990.reproject(dst_res=60)
print(dem_test.info())

# #### **Question:** What is the new raster height?

# #### Change extent/bounds

dem_test = dem_1990.reproject(dst_bounds={"left":502810, "top":8674000, "right":529430, "bottom": 8654290})
print(dem_test.info())

# #### **Question:** What is the new raster height?

# #### Change Coordinate Reference System (CRS) i.e. projection

dem_test = dem_1990.reproject(dst_crs='epsg:4326')
print(dem_test.info())

# ### Reproject the outlines in the same coordinate system as DEMs
# #### **Question:** What is the new pixel resolution? What are the units?

outlines_proj = gu.Vector(outlines_1990.ds.to_crs(dem_2009.crs))
# ## Reproject the outlines in the same coordinate system as DEMs

# ### Calculate the elevation change
outlines_proj = outlines_1990.reproject(dst_crs=dem_2009.crs)

# # Calculating the difference between two DEMs

ddem = dem_2009 - dem_1990

# ### Plot the elevation change map
# ## Plot the elevation change map
# #### Note:
# - `ax` is used here the share the same subplot between the raster and outlines (the default is to create a new figure)
# - `ax` is used here to share the same subplot between the raster and outlines (the default is to create a new figure)
# - ddem is plotted last, to preserve the extent, as glacier outlines cover all of Svalbard.
# - zorder is used to plot in the right sequence (outlines on top)
# - we save the figure to a png file with command `plt.savefig`. You can use a different dpi setting to change the image resolution/size.

# +
vmax = max(abs(np.max(ddem.data)), abs(np.min(ddem.data)))
# Here we automatically calculate the min and max value for the color scale (vmin, vmax). We take the maximum absolute elevation change value. Since the colorbar is divergent (red for negative values, blue for positive) it only makes sense if centered on 0, so we set vmin=-vmax.

ax = plt.subplot(111)
outlines_proj.ds.plot(ax=ax, facecolor='none', edgecolor='k', zorder=2)
ddem.show(ax=ax, cmap='RdYlBu', vmin=-vmax, vmax=vmax, cb_title='Elevation change 2009 - 1990 (m)', zorder=1)
vmax = np.max(np.abs(ddem.data))

fig, ax = plt.subplots(figsize=(10, 8))
outlines_proj.show(ax=ax, facecolor='none', edgecolor='k', zorder=2)
ddem.show(ax=ax, cmap='RdYlBu', vmin=-vmax, vmax=vmax, cbar_title='Elevation change 2009 - 1990 (m)', zorder=1)
ax.set_title('Thinning glaciers near Longyearbyen')
plt.tight_layout()
plt.savefig("ddem_map.png")
plt.show()

# #### <span style='color:red '> **TO DO:** </span> Make changes to the plot:
# 1) Plot the glacier outlines in red.
# 2) Use a different min/max value for colorscale.
# 3) Use a different colormap (choose among the list [here](https://matplotlib.org/stable/gallery/color/colormap_reference.html)).
# 4) Change the figure title.
#
# (Uncomment the lines below and replace "..." with a different value)

# +
# fig, ax = plt.subplots(figsize=(10, 8))
# outlines_proj.show(ax=ax, facecolor='none', edgecolor='...', zorder=2)
# ddem.show(ax=ax, cmap='coolwarm', vmin=..., vmax=..., cbar_title='Elevation change 2009 - 1990 (m)', zorder=1)
# ax.set_title('...')
# plt.tight_layout()
# plt.show()
# -

# ### Saving the results
# ## Saving the elevation change to GTiff

# If you want to export the results to a GTiff for archiving or analyzing in another software (e.g. QGIS), here's how to do it.

ddem.save("temp_ddem.tif")

# ### Rasterize the glacier outlines on the same grid as ddem
# # Calculate zonal statistics
# Here we want to calculate the mean elevation change on and off glaciers.

# ## Rasterize the glacier outlines on the same grid as ddem
# `glacier_mask` is `True` on glaciers, `False` elsewhere.

glacier_mask = outlines_1990.create_mask(ddem)
plt.imshow(glacier_mask.squeeze(), interpolation='none')
glacier_mask.show(add_cbar=False)

# ### Calculate mean dh over glaciers or stable terrain

# Over glaciers:

print(np.mean(ddem.data[glacier_mask]))
print(np.mean(ddem[glacier_mask]))

# Over stable terrain

print(np.mean(ddem.data[~glacier_mask]))
print(np.mean(ddem[~glacier_mask]))



# ### Something is incorrect... mean dh over stable terrain should be ~0 => we need to coregister the DEMs. This is some proper work for our next example !


# Something is wrong, mean dh over stable terrain should be ~0 => we need to coregister the DEMs.
# This is some proper work for our next example !

0 comments on commit 3290d1d

Please sign in to comment.