Introduction

In today’s lab, we will be doing our work inside a web browser (if you are connecting from home you will still need the VPN); using a website https://jupyter.gis.unbc.ca.

This is the first of several labs this semester where we will be working primarily with code as opposed to a Graphical User Interface (GUI). There are many advantages to this approach including but not limited to ease of automation, scale-out-processing, ability to dynamically produce results, and being able to work close to the data.

  • Automation: While it is true that many GIS / Remote Sensing software do include some form of graphical model builder, many also provide this automation through Aplication Interfaces (APIs), with python being a particularly common choice, for example ArcGIS has ArcPy, Qgis QPy, and GDAL provides Python bindings.
  • Scale-out-processing: When working with large datasets that can be divided into many jobs (ie calculating NDVI on 100 images could be 100 jobs) we can use code to coordinate distributing these jobs to sets of computers.
  • Dynamic Results: If we are able to produce a workflow for generating a product, perhaps even through experimentation in a GUI such as PCI Catalys, if we are able to convert these steps to code, we can configure a server to repeat these operations as new data becomes available, or process speific scenes based on user input.
  • Ability to work close to data: While we absolutly could run GUI software anywhere in the world through remote desktop, is it more typical to use code for this. And the motivation is simple, if you need NDVI, that is a relativly small dataset compared to a Sentinal scene. Rather than download the entire scene and then produce NDVI, what if we could produce NDVI in the datacenter where the imagery is stored and download only the NDVI? (Take a look at why and how Digital Globe moved their data to Amazon AWS: https://aws.amazon.com/solutions/case-studies/digitalglobe/)

What is Jupyter

In today’s lab we are using an opensource software stack called JupyterHub and it is made of several components

JupyterHub: This is the webpage you are greeted with when going to https://jupyter.gis.unbc.ca. Once you log in it will spawn a JupyterLab for you to use. In some installations, it may provide different labs to different users, where there is different software and or hardware available. For example the Microsoft Planetary Computer we signed up for in Lab 1 is a JupyterHub, and provides the following JupyterLabs

JupyterLab: This is a personal virtual environment just for you! In our case, we are all using the jupyter-scipy-notebook with some extra geospatial stuff installed. And rather than a personal allotment of CPU and RAM we are all sharing 24 cores and 128GiB of memory.

https://jupyter-docker-stacks.readthedocs.io/en/latest/using/selecting.html#jupyter-scipy-notebook

JupyterNotebook: This is the code editor that we are using, and allows us to write and run our code in the web browser. We can install JupyterNotebook directly on our own computers if we want to (an easy way is with Anaconda Personal: https://www.anaconda.com/products/individual). If we were not in the web browser we could also use Integrated Development Environments such as PyCharm or Spyder. Or text editors such as VSCode, Vim, or Notepad then run our code on the command line.

Anaconda: Anaconda is the programming language that we are using today, it is a flavour of Python.

Lets Get Started

Important note! The JupyterLab environment is temporary and files not moved to the gis file system will be deleted!!!

Log into: https://jupyter.gis.unbc.ca

Here you are presented with the JupyterLab Homescreen. Before we begin we will need access to the GIS Filesystem, under Other select the Terminal and you should see a window as below

The initial text we see tells us the active Anaconda Environment (out of the scope of this class but if you are curious: https://docs.conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html), next is the username of the logged in user…. but who is jovyan???? and finally @ the computer name, notice how this is just a random string as if this computer isn’t really configured for us personally.

enter the command

gisfs/mount.sh

It will then ask for your UNBC username and password enter these there should be no output if the mount is successful, if it says permission denied you likely had a typo in your password, try again and if it still does not work ask your TA for assistance.

To test enter the following command

ls gisfs/K

Do you see a list of your files? If not ask your TA for help.

Now that your files are mounted let’s make a notebook.

To get back to the home screen go to File >>> New Launcher or press [Ctrl] + [Shift] + [L]

Using the file browser in the left side pane, navigate into your GEOG457 folder on K (preferably in a lab3 subfolder), then choose Python 3 (ipykernel) under Notebook. Note this will make Untitiled.ipynb in the active folder to rename simply right click and choose rename from the options.

Coregistration

We will be using the libraries Xdem to co-regiser the DEM’s and then matplotlib to display our results.

In the first block enter the following code

print("Import Libraries...")
import matplotlib.pyplot as plt
import xdem
import numpy as np
import geoutils as gu
print("Done")

Then press the play button

While the code is running the [1] will show as [*] and the results will be displayed below.

Now let us make it easier to find our files, in the next cell make a variable to where the data is stored

src_dir = '/home/jovyan/gisfs/L/GEOG457/Lab3/'
print("Source Directory is " + src_dir)

Then we will open the first DEM, and put it into a variable called reference_dem

print("Load Reference Data...")
reference_dem = xdem.DEM(src_dir + '20170928_columbia_5m_subset.tif')
print("Reference Coordinate Reference: " + str(reference_dem.crs))
print("Nodata vale: " + str(reference_dem.nodata))
reference_dem.data[reference_dem.data==reference_dem.nodata]=np.nan
print("Done")

Before we can plot the image we will calculate the image extents

plt_extent = [
    reference_dem.bounds.left,
    reference_dem.bounds.right,
    reference_dem.bounds.bottom,
    reference_dem.bounds.top,
]

And take a look at the DEM

plt.figure(figsize=(15, 12))
plt.imshow(reference_dem.data.squeeze(), extent=plt_extent)
plt.colorbar()
plt.show()

As we are looking at Glaciers next we will open a shapefile outlining the glaciers, and create a raster mask, then visualize it.

glacier_outlines = gu.Vector(src_dir + 'ice2.shp')
print("Outline CRS: " + str(glacier_outlines.crs))
inlier_mask = ~glacier_outlines.create_mask(reference_dem)
plt.figure(figsize=(15, 12))
plt.imshow(inlier_mask.squeeze(), extent=plt_extent)
plt.show()

The final file we need to open is the Unregisered DEM

dem_to_be_aligned = xdem.DEM(src_dir + '20200828_columbia_5m_subset.tif').reproject(reference_dem, silent=True, resampling='bilinear')
dem_to_be_aligned.data[dem_to_be_aligned.data==dem_to_be_aligned.nodata]=np.nan
plt.figure(figsize=(15, 12))
plt.imshow(dem_to_be_aligned.data.squeeze(), extent=plt_extent)
plt.show()

Before Coregistration lets take a look at the difference between these DEM’s

diff_before = (dem_to_be_aligned - reference_dem).data

plt.figure(figsize=(15, 12))
plt.imshow(diff_before.squeeze(), cmap="coolwarm_r", vmin=-10, vmax=10, extent=plt_extent)
plt.colorbar()
plt.show()

And now for the actual work apply coregistration, note that this may take a minute to process

nuth_kaab = xdem.coreg.NuthKaab()
nuth_kaab.fit(reference_dem.data, dem_to_be_aligned.data, inlier_mask, transform=reference_dem.transform)
aligned_dem_data = nuth_kaab.apply(dem_to_be_aligned.data, transform=dem_to_be_aligned.transform)

Now lets take a look at the change after coregistration

diff_after = aligned_dem_data - reference_dem.data

plt.figure(figsize=(15, 12))
plt.imshow(diff_after.squeeze(), cmap="coolwarm_r", vmin=-10, vmax=10, extent=plt_extent)
plt.colorbar()
plt.show()
error_before = diff_before + 0 #Plus 0 prevents damage to diff_before
error_before[~inlier_mask] = np.nan

error_after = diff_after + 0 #Plus 0 prevents damage to diff_before
error_after[~inlier_mask] = np.nan

print(f"Error before: {xdem.spatialstats.nmad(error_before):.2f} m")
print(f"Error after: {xdem.spatialstats.nmad(error_after):.2f} m")

And just for fun, the difference between the resutls

diff_diff = diff_before - diff_after

plt.figure(figsize=(15, 12))
plt.imshow(diff_diff.squeeze(), cmap="coolwarm_r", vmin=-3, vmax=3, extent=plt_extent)
plt.colorbar()
plt.show()

Assignment

This weeks assignment is worth 5 points

  • (0.5) Why should we coregister DEM’s before calculating changes?
  • (0.5) Why do we use a Glacier Mask when coregistering these DEM’s?
  • (0.5) You do not need to implement, but given a scenario where you are not provided with a glacier mask, how might you produce a mask using Remote Sensing? Assume you have access to multispectral data, and that you do not wish to trace the glacier by hand. (expectation of detail is 1-2 sentences)
  • (0.5) If we wanted to look at changes in other environments, what are some things we may wish to mask? Provide at least one example in an urban setting and one in a non-urban setting
  • Perform coregistration on DEM’s located in L:\GEOG457\Lab3\monkman; use trim1985clip3.tif as the reference alos2005float.tif as to be aligned, and glaciers2013.shp for the mask. Note there are lots of files here. Make sure you are using the ones specified.
    • (0.5) Ensure all plots have adequate contrast before submitting
    • (0.5) Change the colours of plots from what was used in the lab; the colours chosen should still be intuitive for viewers of your plot. (There is great documentation on Matplotlib.org)
    • (0.5) In Science, communication is essential; what guidance does Matlibplot provide for choosing colour schemes accessible to those with colour vision deficiencies?
    • (0.5) Calculate the change in volume from 1985 to 2005. As a hint, see bottom of page for code that would work for the DEM’s that we used in Lab; however, there is a catch the data we are using in the assignment has a different format, while the X, Y is still in metres, the Z values are not; even in Canada older DEMs use feet for elevation. Modify the code to give a correct volume change.
    • (1) Please submit your completed notebook (File >>> Download), put it into a zip folder before sending it as UNBC’s email system will block the .ipynb file extension
pixel_area = reference_dem.res[0] * reference_dem.res[1] # X and Y values of resolution
#volume = sum of elevation not in inlier_mask times pixel area
#1000000000m^3 in 1km^3
nonreg_loss = (np.nansum(diff_before[~inlier_mask]) * pixel_area) / 1000000000
print("Change in Glacier volume without coregistration: " + str(nonreg_loss) + "km^3")
reg_loss = (np.nansum(diff_after[~inlier_mask]) * pixel_area) / 1000000000
print("Change in Glacier volume with coregistration: " + str(reg_loss) + "km^3")
Categories: GEOG 457Labs