.ipynb .pdf

repository open issue suggest edit

Binder

Contents

Outputs export

Several utilities are available for exporting model outputs in useful formats.

We will go through the main ones.

import os
import pooch
from script import readOutput as rout

import matplotlib
import pyvista as pv
import matplotlib.pyplot as plt

label_size = 7
matplotlib.rcParams['xtick.labelsize'] = label_size
matplotlib.rcParams['ytick.labelsize'] = label_size
matplotlib.rc('font', size=6)

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

We first define a folder where exported files will be stored:

out_path = 'export'
if not os.path.exists(out_path):
    os.makedirs(out_path)

Netcdf outputs

We start with the netcdf outputs as it is the most common one used in our field.

Building netcdf file

Netcdf exports are done by using the readOutput.py script presented in the previous notebook. Here we export all the time steps at once by looping through the number of outputs (5 in this case).

Arguments for readOutput.py

The readOutput.py script main function requires several arguments:

• path: the path to the input file

• filename: the name of the input file

• step: the step you wish to output (here set to 5 corresponding to the last output based on the input parameters: start time 0 year, end time 50 thousand years with an output every 10 thousand years)

• nbstep: the number of time steps to plot (useful if one want to output a netdcf file containing all time steps (done in the following section).

• uplift_forcing: set to False as we are not considering any tectonic forcing

Then the buildLonLatMesh function is used to interpolate (using a kd-tree approach) the gospl variables on a regular mesh. It also provides a way to limit the created netcdf file by defining a bounding box:

uplift_forcing = False

# Specifying the grid resolution in degrees
reso = 0.1

# Total number of outputs
nbstep = 5

# Bounding box
bb = [-134,20,-46,80]

def build_ncgrid(ncout, nbstep, bbox, reso, fname, uplift_forcing):

    # Looping through the output time steps
    for k in range(nbstep+1):
        if k == 0:
            # Calling the initialisation function for our class 
            ncgrid = rout.readOutput(path='./', filename=fname, 
                                     step=k, nbstep=nbstep+1, 
                                     uplift=uplift_forcing)
        else:
            # Update the variables after the first time steps
            ncgrid.update(step=k, uplift=uplift_forcing)

        # Build the regular grid defining the bounding box
        ncgrid.buildLonLatMesh(res=reso, nghb=3, box=bbox)
    
    # Export as netcdf
    ncgrid.exportNetCDF(ncfile = ncout)
    
    return

Calling the function build_ncgrid to export the netcdf file on the desired bounding box:

ncout = os.path.join(out_path, "GoMresult.nc")
# We commented the next line and loaded the dataset from figshare (see below)
#build_ncgrid(ncout, nbstep, bb, reso, fname='input.yml', uplift_forcing=False)

Here we will download the file directly from figshare:

ncout = pooch.retrieve(
    url="https://ndownloader.figshare.com/files/27579284",
    known_hash="sha256:fdc4a4ac70bca66ed20c3a126fdbf8fd2b33a2b47170e45c53eaf96f20f42a32",
    downloader=pooch.HTTPDownloader(progressbar=True),
    fname="GoMresult.nc",
)

Downloading data from 'https://ndownloader.figshare.com/files/27579284' to file '/home/runner/.cache/pooch/GoMresult.nc'.

  0%|                                              | 0.00/74.8M [00:00<?, ?B/s]

  0%|                                      | 52.2k/74.8M [00:00<03:15, 383kB/s]

  0%|                                       | 157k/74.8M [00:00<02:05, 593kB/s]

  1%|▎                                     | 557k/74.8M [00:00<00:45, 1.62MB/s]

  3%|█                                    | 2.21M/74.8M [00:00<00:12, 5.60MB/s]

  9%|███▎                                 | 6.70M/74.8M [00:00<00:04, 15.0MB/s]

 15%|█████▍                               | 11.0M/74.8M [00:00<00:03, 20.3MB/s]

 21%|███████▋                             | 15.6M/74.8M [00:00<00:02, 26.8MB/s]

 25%|█████████▍                           | 19.1M/74.8M [00:01<00:01, 27.9MB/s]

 31%|███████████▌                         | 23.4M/74.8M [00:01<00:01, 32.1MB/s]

 36%|█████████████▎                       | 26.9M/74.8M [00:01<00:01, 31.7MB/s]

 42%|███████████████▌                     | 31.4M/74.8M [00:01<00:01, 35.3MB/s]

 47%|█████████████████▎                   | 35.1M/74.8M [00:01<00:01, 33.9MB/s]

 53%|███████████████████▌                 | 39.6M/74.8M [00:01<00:00, 37.1MB/s]

 58%|█████████████████████▍               | 43.4M/74.8M [00:01<00:00, 35.4MB/s]

 65%|███████████████████████▉             | 48.3M/74.8M [00:01<00:00, 36.1MB/s]

 71%|██████████████████████████▎          | 53.3M/74.8M [00:01<00:00, 39.6MB/s]

 77%|████████████████████████████▎        | 57.3M/74.8M [00:02<00:00, 37.7MB/s]

 83%|██████████████████████████████▊      | 62.3M/74.8M [00:02<00:00, 40.9MB/s]

 89%|████████████████████████████████▊    | 66.4M/74.8M [00:02<00:00, 38.3MB/s]

 94%|██████████████████████████████████▊  | 70.5M/74.8M [00:02<00:00, 36.7MB/s]

  0%|                                              | 0.00/74.8M [00:00<?, ?B/s]

100%|█████████████████████████████████████| 74.8M/74.8M [00:00<00:00, 75.0GB/s]

Visualisation with ipygany

Visualise the result in Jupyter environment with ipygany. This is done using by opening the netcdf file first:

import numpy as np
import xarray as xr

from ipygany import PolyMesh, Scene, IsoColor, WarpByScalar, ColorBar, colormaps
from ipywidgets import VBox, FloatSlider, FileUpload, Dropdown, jslink

ds = xr.open_dataset(ncout, decode_times=False)
ds

<xarray.Dataset>
Dimensions:        (latitude: 681, longitude: 881, time: 6)
Coordinates:
  * time           (time) float64 0.0 1e+04 2e+04 3e+04 4e+04 5e+04
  * latitude       (latitude) float64 12.0 12.1 12.2 12.3 ... 79.8 79.9 80.0
  * longitude      (longitude) float64 -134.0 -133.9 -133.8 ... -46.1 -46.0
Data variables:
    elevation      (time, latitude, longitude) float64 ...
    erodep         (time, latitude, longitude) float64 ...
    precipitation  (time, latitude, longitude) float64 ...
    drainageArea   (time, latitude, longitude) float64 ...
    basinID        (time, latitude, longitude) int32 ...
Attributes:
    description:  gospl outputs
    history:      Created Thu Apr  8 17:39:07 2021

We then select a specific time step and variable using xarray functions:

# Selecting last time step
ds_z = ds.isel(time=[-1])

# Dropping all variables expect the elevation
ds_z = ds_z.drop(['time','erodep', 'precipitation', 'drainageArea', 'basinID'])
ds_z = ds_z.squeeze("time")
ds_z

<xarray.Dataset>
Dimensions:    (latitude: 681, longitude: 881)
Coordinates:
  * latitude   (latitude) float64 12.0 12.1 12.2 12.3 ... 79.7 79.8 79.9 80.0
  * longitude  (longitude) float64 -134.0 -133.9 -133.8 ... -46.2 -46.1 -46.0
Data variables:
    elevation  (latitude, longitude) float64 ...
Attributes:
    description:  gospl outputs
    history:      Created Thu Apr  8 17:39:07 2021

Note

We now create pyvista structured mesh (our netcdf is structured!)

xx, yy, zz = np.meshgrid(np.radians(ds_z['longitude']), 
                         np.radians(ds_z['latitude']), 
                         [0])

# Transform to spherical coordinates
radius = 6371.0e6
x = radius * np.cos(yy) * np.cos(xx)
y = radius * np.cos(yy) * np.sin(xx)
z = radius * np.sin(yy)

grid = pv.StructuredGrid(x, y, z)

# Add data to mesh
for var in ds_z.data_vars:
    grid[var] = np.array(ds_z[var]).ravel(order='F')

We can then convert pyvista mesh to ipygany mesh

See also

This approach is based on the following cerege example.

# Turn the PyVista mesh into a PolyMesh
mesh = PolyMesh.from_vtk(grid)

# Color the mesh
colored_mesh = IsoColor(mesh, min=ds_z.elevation.min(), max=ds_z.elevation.max())

# setup warping
warped_mesh = WarpByScalar(colored_mesh, input='elevation', factor=0)

# Link a slider to the warp value
warp_slider = FloatSlider(min=0., max=10., value=5)

def on_slider_change(change):
    warped_mesh.factor = change['new'] * -10000

warp_slider.observe(on_slider_change, 'value')

# Create a colorbar widget
colorbar = ColorBar(colored_mesh)

# Colormap choice widget
colormap = Dropdown(
    options=colormaps,
    description='colormap:'
)

jslink((colored_mesh, 'colormap'), (colormap, 'index'))

VBox((colormap, warp_slider, Scene([warped_mesh])))

Geotiff output

xarray functionality

Let first use the netcdf file created and open it with the xarray library in the jupyter environment:

import xarray as xr

ds = xr.open_dataset(ncout, decode_times=False)
ds

<xarray.Dataset>
Dimensions:        (latitude: 681, longitude: 881, time: 6)
Coordinates:
  * time           (time) float64 0.0 1e+04 2e+04 3e+04 4e+04 5e+04
  * latitude       (latitude) float64 12.0 12.1 12.2 12.3 ... 79.8 79.9 80.0
  * longitude      (longitude) float64 -134.0 -133.9 -133.8 ... -46.1 -46.0
Data variables:
    elevation      (time, latitude, longitude) float64 ...
    erodep         (time, latitude, longitude) float64 ...
    precipitation  (time, latitude, longitude) float64 ...
    drainageArea   (time, latitude, longitude) float64 ...
    basinID        (time, latitude, longitude) int32 ...
Attributes:
    description:  gospl outputs
    history:      Created Thu Apr  8 17:39:07 2021

By default the mesh is written in lon/lat (projection epsg:4326 as gospl is a global model).

Using the rioxarray library we have the ability to reproject the dataset in any other type of projection. Let’s reproject the dataset in utm coordinates:

import rioxarray

ds = ds.rio.write_crs(4326)
print('Default projection:',ds.rio.crs)

print('Estimated UTM projection:',ds.rio.estimate_utm_crs())
ds_utm = ds.rio.reproject(ds.rio.estimate_utm_crs())

ds_utm

Default projection: EPSG:4326
Estimated UTM projection: EPSG:32615

<xarray.Dataset>
Dimensions:        (time: 6, x: 1255, y: 937)
Coordinates:
  * x              (x) float64 -4.36e+06 -4.352e+06 ... 6.225e+06 6.233e+06
  * y              (y) float64 9.233e+06 9.224e+06 ... 1.334e+06 1.326e+06
  * time           (time) float64 0.0 1e+04 2e+04 3e+04 4e+04 5e+04
    spatial_ref    int64 0
Data variables:
    elevation      (time, y, x) float64 1.798e+308 1.798e+308 ... 1.798e+308
    erodep         (time, y, x) float64 1.798e+308 1.798e+308 ... 1.798e+308
    precipitation  (time, y, x) float64 1.798e+308 1.798e+308 ... 1.798e+308
    drainageArea   (time, y, x) float64 1.798e+308 1.798e+308 ... 1.798e+308
    basinID        (time, y, x) int32 -2147483647 -2147483647 ... -2147483647
Attributes:
    description:  gospl outputs
    history:      Created Thu Apr  8 17:39:07 2021

Let’s now create a geotiff file containing the elevation for the last time step using the rioxarray functionality:

elevArray = ds_utm.isel(time=[-1])["elevation"][0,:,:]

# Export to geotiff
tifout = os.path.join(out_path, "GoMresult5.tif")

elevArray.rio.to_raster(tifout)

Advanced functionalities on geotiff

We can use the rasterio library to visualise the geotiff file in our notebook:

import rasterio
from rasterio.mask import mask
from rasterio.plot import show, show_hist

data = rasterio.open(tifout)

print('Number of band',data.count)
print('Image resolution',data.height, data.width)
print('CRS',data.crs)

Number of band 1
Image resolution 937 1255
CRS EPSG:32615

fig, (axr) = plt.subplots(1,1, figsize=(8,9))
show((data, 1), ax=axr, cmap='gist_earth', vmin=-10000, vmax=10000)
plt.show()

What we want to do next is to create a bounding box around the Gulf of Mexico region and clip the raster based on that.

We create a bounding box with Shapely.

from shapely.geometry import box

# WGS84 coordinates (lon/lat - x/y)
minx, miny = -102.5, 18.5
maxx, maxy = -82.0, 40.0 
bbox = box(minx, miny, maxx, maxy)

We insert the bbox into a GeoDataFrame and re-project into the same coordinate system as the raster data

import geopandas as gpd
from fiona.crs import from_epsg

geo = gpd.GeoDataFrame({'geometry': bbox}, index=[0], crs=from_epsg(4326))
geo = geo.to_crs(crs=data.crs.data)

Next we need to get the coordinates of the geometry in such a format that rasterio wants them. This can be conducted easily with following function

def getFeatures(gdf):
    
    """Function to parse features from GeoDataFrame in such a manner that rasterio wants them"""
    
    import json
    return [json.loads(gdf.to_json())['features'][0]['geometry']]

Get the geometry coordinates by using the function.

coords = getFeatures(geo)
print(coords)

[{'type': 'Polygon', 'coordinates': [[[1666930.9647581263, 2081361.9678710762], [1439944.6735851327, 4486147.723400535], [-311558.0298494082, 4471222.562471898], [-506532.72480437916, 2072156.005210298], [1666930.9647581263, 2081361.9678710762]]]}]

Now we are ready to clip the raster with the polygon using the coords variable that we just created.

Tip

Clipping the raster can be done easily with the mask function that we imported in the beginning from rasterio, and specifying clip=True.

out_img, out_transform = mask(data, shapes=coords, crop=True)

Next, we need to modify the metadata.

Let’s start by copying the metadata from the original data file.

Then we parse the EPSG value from the CRS so that we can create a Proj4 string using PyCRS library (to ensure that the projection information is saved correctly).

import pycrs

out_meta = data.meta.copy()
print(out_meta)

epsg_code = int(data.crs.data['init'][5:])
print(epsg_code)


out_meta.update({"driver": "GTiff", "height": out_img.shape[1],
                 "width": out_img.shape[2], "transform": out_transform,
                 "crs": pycrs.parse.from_epsg_code(epsg_code).to_proj4()})

{'driver': 'GTiff', 'dtype': 'float64', 'nodata': 1.7976931348623157e+308, 'width': 1255, 'height': 937, 'count': 1, 'crs': CRS.from_epsg(32615), 'transform': Affine(8447.814056310064, 0.0, -4364578.0519712195,
       0.0, -8447.814056310064, 9237021.067040086)}
32615

Finally, we can save the clipped raster to disk with following command.

tifout2 = os.path.join(out_path, "GoM_clipped.tif")

with rasterio.open(tifout2, "w", **out_meta) as dest:
    dest.write(out_img)

Let’s still check that the result is correct by plotting our new clipped raster.

clipped = rasterio.open(tifout2)

fig, (axr) = plt.subplots(1,1, figsize=(5,5))
show((clipped, 1), ax=axr, cmap='gist_earth', vmin=-1000, vmax=1000)
plt.show()

ZMAP files

The zmapio library allows to read and write map gridded data using ZMAP Plus ASCII Grid format.

Here we will use it to write our elevation as a ZMAP grid.

from zmapio import ZMAPGrid

z_val = elevArray.values[::5,::5].T
print('Z-values shape: ', z_val.shape)

Z-values shape:  (251, 188)

Define the ZMAP dataset:

new_zgrid = ZMAPGrid(z_values=z_val, min_x=ds_utm.x.values.min(), 
                     max_x=ds_utm.x.values.max(),
                     min_y=ds_utm.y.values.min(),  
                     max_y=ds_utm.y.values.max())

Write the new ZMAP file by customising the formating:

zgridout = os.path.join(out_path, "GoM_zmap.dat")

new_zgrid.comments = ['Model', 'output']
new_zgrid.nodes_per_line = 4
new_zgrid.field_width = 15
new_zgrid.decimal_places = 3
new_zgrid.name = 'gospl'
new_zgrid.write(zgridout)

!head $zgridout

!Model
!output
@gospl, GRID, 4
15, 1e+30, , 3, 1
188, 251, -4360354.144943064, 6233204.681669757, 1325643.203305711, 9232797.16001193
0.0, 0.0, 0.0
@
179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000
179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000
179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000

Let’s visualise the dataset…

new_zgrid.plot(cmap='gist_earth', shading='auto')

<matplotlib.collections.QuadMesh at 0x7f661e6d2b20>

Landscape evolution Catchment characteristics