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 filefilename
: the name of the input filestep
: 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 anetdcf
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
- latitude: 681
- longitude: 881
- time: 6
- time(time)float640.0 1e+04 2e+04 3e+04 4e+04 5e+04
- units :
- days since 0000-01-01
- calendar :
- 365_day
array([ 0., 10000., 20000., 30000., 40000., 50000.])
- latitude(latitude)float6412.0 12.1 12.2 ... 79.8 79.9 80.0
- units :
- degrees_north
array([12. , 12.1, 12.2, ..., 79.8, 79.9, 80. ])
- longitude(longitude)float64-134.0 -133.9 ... -46.1 -46.0
- units :
- degrees_east
array([-134. , -133.9, -133.8, ..., -46.2, -46.1, -46. ])
- elevation(time, latitude, longitude)float64...
- units :
- metres
[3599766 values with dtype=float64]
- erodep(time, latitude, longitude)float64...
- units :
- metres
[3599766 values with dtype=float64]
- precipitation(time, latitude, longitude)float64...
- units :
- m/yr
[3599766 values with dtype=float64]
- drainageArea(time, latitude, longitude)float64...
- units :
- m2
[3599766 values with dtype=float64]
- basinID(time, latitude, longitude)int32...
- units :
- int
[3599766 values with dtype=int32]
- 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
- latitude: 681
- longitude: 881
- latitude(latitude)float6412.0 12.1 12.2 ... 79.8 79.9 80.0
- units :
- degrees_north
array([12. , 12.1, 12.2, ..., 79.8, 79.9, 80. ])
- longitude(longitude)float64-134.0 -133.9 ... -46.1 -46.0
- units :
- degrees_east
array([-134. , -133.9, -133.8, ..., -46.2, -46.1, -46. ])
- elevation(latitude, longitude)float64...
- units :
- metres
[599961 values with dtype=float64]
- 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
# 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
- latitude: 681
- longitude: 881
- time: 6
- time(time)float640.0 1e+04 2e+04 3e+04 4e+04 5e+04
- units :
- days since 0000-01-01
- calendar :
- 365_day
array([ 0., 10000., 20000., 30000., 40000., 50000.])
- latitude(latitude)float6412.0 12.1 12.2 ... 79.8 79.9 80.0
- units :
- degrees_north
array([12. , 12.1, 12.2, ..., 79.8, 79.9, 80. ])
- longitude(longitude)float64-134.0 -133.9 ... -46.1 -46.0
- units :
- degrees_east
array([-134. , -133.9, -133.8, ..., -46.2, -46.1, -46. ])
- elevation(time, latitude, longitude)float64...
- units :
- metres
[3599766 values with dtype=float64]
- erodep(time, latitude, longitude)float64...
- units :
- metres
[3599766 values with dtype=float64]
- precipitation(time, latitude, longitude)float64...
- units :
- m/yr
[3599766 values with dtype=float64]
- drainageArea(time, latitude, longitude)float64...
- units :
- m2
[3599766 values with dtype=float64]
- basinID(time, latitude, longitude)int32...
- units :
- int
[3599766 values with dtype=int32]
- 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
- time: 6
- x: 1255
- y: 937
- x(x)float64-4.36e+06 -4.352e+06 ... 6.233e+06
- axis :
- X
- long_name :
- x coordinate of projection
- standard_name :
- projection_x_coordinate
- units :
- metre
array([-4360354.144943, -4351906.330887, -4343458.51683 , ..., 6216309.053557, 6224756.867613, 6233204.68167 ])
- y(y)float649.233e+06 9.224e+06 ... 1.326e+06
- axis :
- Y
- long_name :
- y coordinate of projection
- standard_name :
- projection_y_coordinate
- units :
- metre
array([9232797.160012, 9224349.345956, 9215901.531899, ..., 1342538.831418, 1334091.017362, 1325643.203306])
- time(time)float640.0 1e+04 2e+04 3e+04 4e+04 5e+04
- units :
- days since 0000-01-01
- calendar :
- 365_day
array([ 0., 10000., 20000., 30000., 40000., 50000.])
- spatial_ref()int640
- crs_wkt :
- PROJCS["WGS 84 / UTM zone 15N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-93],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32615"]]
- semi_major_axis :
- 6378137.0
- semi_minor_axis :
- 6356752.314245179
- inverse_flattening :
- 298.257223563
- reference_ellipsoid_name :
- WGS 84
- longitude_of_prime_meridian :
- 0.0
- prime_meridian_name :
- Greenwich
- geographic_crs_name :
- WGS 84
- horizontal_datum_name :
- World Geodetic System 1984
- projected_crs_name :
- WGS 84 / UTM zone 15N
- grid_mapping_name :
- transverse_mercator
- latitude_of_projection_origin :
- 0.0
- longitude_of_central_meridian :
- -93.0
- false_easting :
- 500000.0
- false_northing :
- 0.0
- scale_factor_at_central_meridian :
- 0.9996
- spatial_ref :
- PROJCS["WGS 84 / UTM zone 15N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-93],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32615"]]
- GeoTransform :
- -4364578.0519712195 8447.814056310064 0.0 9237021.067040086 0.0 -8447.814056310064
array(0)
- elevation(time, y, x)float641.798e+308 ... 1.798e+308
- units :
- metres
- _FillValue :
- 1.7976931348623157e+308
array([[[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], ..., [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]], [[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], ... [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]], [[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], ..., [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]]])
- erodep(time, y, x)float641.798e+308 ... 1.798e+308
- units :
- metres
- _FillValue :
- 1.7976931348623157e+308
array([[[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], ..., [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]], [[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], ... [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]], [[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], ..., [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]]])
- precipitation(time, y, x)float641.798e+308 ... 1.798e+308
- units :
- m/yr
- _FillValue :
- 1.7976931348623157e+308
array([[[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], ..., [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]], [[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], ... [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]], [[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], ..., [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]]])
- drainageArea(time, y, x)float641.798e+308 ... 1.798e+308
- units :
- m2
- _FillValue :
- 1.7976931348623157e+308
array([[[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], ..., [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]], [[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], ... [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]], [[1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], ..., [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308], [1.79769313e+308, 1.79769313e+308, 1.79769313e+308, ..., 1.79769313e+308, 1.79769313e+308, 1.79769313e+308]]])
- basinID(time, y, x)int32-2147483647 ... -2147483647
- units :
- int
- _FillValue :
- -2147483647
array([[[-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], ..., [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647]], [[-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], ... [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647]], [[-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], ..., [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647], [-2147483647, -2147483647, -2147483647, ..., -2147483647, -2147483647, -2147483647]]], dtype=int32)
- 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>