Load Sample Source

Illustrate how to load a netCDF file for one source using xarray or rasterio, and how to grab the last frame to make an "instant" rupture version.

Also write it out as a GeoClaw dtopo file with ASCII raster format, as described in the GeoClaw Documentation, and make some plots using GeoClaw tools.

Adapt this code to write it out to whatever format you need.

%matplotlib inline

from pylab import *
from pathlib import Path

dtopo_dir = './dtopofiles_nc'  # path to unzipped directory

event = 'BL13M'
path = Path(dtopo_dir) / f'{event}.nc'
print(path)

dtopofiles_nc/BL13M.nc

Open with xarray:

Using xarray.

Also illustrate how to grab the final static deformation.

import xarray as xr
with xr.open_dataset(path, decode_timedelta=False) as dtopo_xr:
    print(dtopo_xr)

    # To grab only the final deformation for use as an "instant" event:
    
    print(f"\nThe dz array has shape {dtopo_xr.variables['dz'].shape}")
    dz_static = dtopo_xr.variables['dz'][-1,:,:]
    print(f"The dz_static array contains the final (static) deformation, with shape {dz_static.shape}")
    print(f"The maximum static uplift is {dz_static.max():.2f} meters")

<xarray.Dataset> Size: 107MB
Dimensions:  (time: 31, lat: 1201, lon: 721)
Coordinates:
  * time     (time) float64 248B 0.0 10.0 20.0 30.0 ... 270.0 280.0 290.0 300.0
  * lat      (lat) float64 10kB 40.0 40.01 40.02 40.02 ... 49.98 49.99 50.0
  * lon      (lon) float64 6kB -128.5 -128.5 -128.5 ... -122.5 -122.5 -122.5
Data variables:
    dz       (time, lat, lon) float32 107MB ...

The dz array has shape (31, 1201, 721)
The dz_static array contains the final (static) deformation, with shape (1201, 721)
The maximum static uplift is 8.15 meters

Open with rasterio:

Or use rasterio, another standard Python package that supports reading netCDF easily.

import rasterio
with rasterio.open(path) as src:
    print(f"Coordinate Reference System (CRS): {src.crs}")
    print(f"Bounds: {src.bounds}")
    print(f"Number of bands: {src.count}")
    print(f"Width/Height: {src.width}x{src.height}")
    
    meta = src.meta
    print('\nsrc.meta = \n')
    for k in meta.keys():
        print(f'{k:<30}{meta[k]}')
    
    tags = src.tags()
    print('\nsrc.tags = \n')
    for k in tags.keys():
        print(f'{k:<30}:  {tags[k]}')

Coordinate Reference System (CRS): None
Bounds: BoundingBox(left=-128.50416666666666, bottom=39.99583333333333, right=-122.4958333333388, top=50.004166666666094)
Number of bands: 31
Width/Height: 721x1201

src.meta = 

driver                        netCDF
dtype                         float32
nodata                        9.969209968386869e+36
width                         721
height                        1201
count                         31
crs                           None
transform                     | 0.01, 0.00,-128.50|
| 0.00,-0.01, 50.00|
| 0.00, 0.00, 1.00|

src.tags = 

dz#long_name                  :  seafloor deformation
dz#units                      :  meters
lat#axis                      :  Y
lat#standard_name             :  latitude
lat#units                     :  degrees_north
lon#axis                      :  X
lon#standard_name             :  longitude
lon#units                     :  degrees_east
NETCDF_DIM_EXTRA              :  {time}
NETCDF_DIM_time_DEF           :  {31,6}
NETCDF_DIM_time_VALUES        :  {0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,190,200,210,220,230,240,250,260,270,280,290,300}
time#axis                     :  T
time#standard_name            :  time
time#units                    :  seconds

Load with GeoClaw

NOTE: Loading netCDF dtopo files will be available in Clawpack v5.15.0 but not in earlier releases.

See GeoClaw Documentation for info on the file formats and other tools available.

from clawpack.geoclaw import dtopotools

try:
    # this only works with the master branch of geoclaw, not yet in any release
    dtopo = dtopotools.DTopography(path, dtopo_type=4)  # 4 ==> netcdf format
except:
    import clawpack
    print(f'Using clawpack version {clawpack.__version__}, direct load of netCDF not supported')
    print('*** The remainder of this notebook will fail!')

Rewrite as GeoClaw ascii file

fname = f'{event}.dtt3'
dtopo.write(fname, dtopo_type=3)

sizeMB = Path(fname).stat().st_size / 1e6
print(f'{fname} has size {sizeMB:.1f} MB')

BL13M.dtt3 has size 176.8 MB

Capture only the final deformation for an "instant" rupture

dz_static = dtopo.dZ
print(f"\nThe dtopo.dZ array has shape {dtopo.dZ.shape}")
dz_static = dz_static = dtopo.dZ[-1,:,:]
print(f"The dz_static array contains the final (static) deformation, with shape {dz_static.shape}")
print(f"The maximum static uplift is {dz_static.max():.2f} meters")

# Create a new dtopo file for an instant rupture at time 1 second:

import copy
dtopo_instant = copy.copy(dtopo)  # so we have the same X and Y arrays
dtopo_instant.dZ = dtopo_instant.dZ[-2:-1, :,:]  # only the last dz
dtopo_instant.times = [1.0]  # desired time of instant rupture
print(f"\nNew dtopo_instant.dZ has shape {dtopo_instant.dZ.shape}" \
        + f" with maximum {dtopo_instant.dZ.max():.1f}")
print(f"Instant displacement with specified times {dtopo_instant.times}")

fname = f'{event}_instant.dtt3'
dtopo_instant.write(fname, dtopo_type=3)
sizeMB = Path(fname).stat().st_size / 1e6
print(f'{fname} has size {sizeMB:.1f} MB')

The dtopo.dZ array has shape (31, 1201, 721)
The dz_static array contains the final (static) deformation, with shape (1201, 721)
The maximum static uplift is 8.15 meters

New dtopo_instant.dZ has shape (1, 1201, 721) with maximum 8.1
Instant displacement with specified times [1.0]
BL13M_instant.dtt3 has size 5.6 MB

Plots of deformation

Using the GeoClaw dtopotools module tools.

coast = load('../topo/CSZ_coast.npy')  # precomputed coastline

# time to plot deformation
#tplot = dtopo.times.max()  # for final static deformation
tplot1 = dtopo.times.max() / 3.
tplot2 = dtopo.times.max()

fig,axs = subplots(1,2,figsize=(10,6))
for ax in axs:
    ax.plot(coast[:,0], coast[:,1], 'g', linewidth=0.9)
    ax.set_aspect(1/cos(45*pi/180))
    ax.set_xlim(-130,-121)
    ax.set_ylim(39,50)

dtopo.plot_dZ_colors(t=tplot1, axes=axs[0], dZ_interval=100, cmax_dZ=10);
axs[0].set_title(f'{event}\nSeafloor deformation dz\nat time t = {tplot1:.1f} seconds');

dtopo.plot_dZ_colors(t=tplot2, axes=axs[1], dZ_interval=100, cmax_dZ=10);
axs[1].set_title(f'{event}\nSeafloor deformation dz\nat time tfinal = {tplot2:.1f} seconds');

To Do: Add an animation.

Vertical deformation at Lagoon Creek

# make interpolating function of (x,y,t):
dtopo_fcn = dtopo.make_function()

# Evaluate at desired location (gauge 1 location, on beach):
xg = -124.102
yg = 41.596
tg = arange(0, dtopo.times.max(), 1)  # for time series covering deformation time, with dt=1 sec
dz = dtopo_fcn(xg, yg, tg)

figure(figsize=(8,3))
plot(tg, dz, 'g')
grid(True)
xlabel('time (seconds)')
ylabel(f'vertical deformation dz (meters)')
title(f'vertical deformation dz at xg={xg:.5f}, yg={yg:.5f}');