In [1]:
import os, sys
import cmocean
import matplotlib.pyplot as plt
import numpy as np
import copy
from numpy import ma # masked arrays
from clawpack.geoclaw import topotools, marching_front
from clawpack.visclaw.plottools import pcolorcells
#from clawpack.visclaw import colormaps
#from clawpack.amrclaw import region_tools
In [2]:
topofile_org = './topo_05.dat'
topofile_new = topofile_org.replace('.dat','_masked.dat')
topofile_new = os.path.basename(topofile_new)
topo = topotools.Topography(topofile_org, topo_type=3)
In [3]:
crange = [-20.0,50.0]

cmap = cmocean.cm.topo
newcmap = cmocean.tools.crop(cmap, crange[0], crange[1], 0)

fig = plt.figure()
ax = fig.add_subplot(111)
pc = pcolorcells(topo.X, topo.Y, topo.Z, cmap=newcmap, vmin=crange[0], vmax=crange[1])
ax.axis('scaled')
cbar = fig.colorbar(pc)
In [4]:
wet_points = marching_front.select_by_flooding(topo.Z, Z1=-4.0, Z2=-0.30, max_iters=None)
Zdry = ma.masked_array(topo.Z, wet_points)

## plot
fig = plt.figure(figsize=(12,8))
# topography with wet points masked
ax1 = fig.add_subplot(121)
pc = ax1.pcolor(topo.X, topo.Y, Zdry, cmap=newcmap, vmin=crange[0], vmax=crange[1])
ax1.axis('scaled')
cbar = fig.colorbar(pc)
# wet points
ax2 = fig.add_subplot(122)
pc = ax2.pcolor(topo.X, topo.Y, wet_points)
ax2.axis('scaled')
cbar = fig.colorbar(pc)

Selecting points with Z1 = -4, Z2 = -0.3, max_iters=586551
Done after 155 iterations with 146286 points chosen

<ipython-input-4-013e1a1d5718>:8: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  pc = ax1.pcolor(topo.X, topo.Y, Zdry, cmap=newcmap, vmin=crange[0], vmax=crange[1])
<ipython-input-4-013e1a1d5718>:13: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  pc = ax2.pcolor(topo.X, topo.Y, wet_points)
In [5]:
## Create a version of topo.Z with all wet points masked out:
mask_dry = np.logical_not(wet_points)
Z_dry = ma.masked_array(topo.Z, wet_points) 

## Create a version of topo.Z with only dry points below MHW masked out:
mask_dry_onshore = np.logical_and(mask_dry, topo.Z<0.)
Z_allow_wet = ma.masked_array(topo.Z, mask_dry_onshore)

## plot
fig = plt.figure(figsize=(12,8))
# Z_allow_wet
ax1 = fig.add_subplot(121)
pc = ax1.pcolor(topo.X, topo.Y, Z_allow_wet, cmap=newcmap, vmin=crange[0], vmax=crange[1])
cbar = fig.colorbar(pc)
ax1.axis('scaled')
# mask_dry_onshore
ax2 = fig.add_subplot(122)
pc = ax2.pcolor(topo.X, topo.Y, mask_dry_onshore, vmin=0, vmax=1)
cbar = fig.colorbar(pc)
ax2.axis('scaled')

<ipython-input-5-355918889f69>:13: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  pc = ax1.pcolor(topo.X, topo.Y, Z_allow_wet, cmap=newcmap, vmin=crange[0], vmax=crange[1])
<ipython-input-5-355918889f69>:18: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  pc = ax2.pcolor(topo.X, topo.Y, mask_dry_onshore, vmin=0, vmax=1)
Out[5]:
(-11730.0, 7770.0, -49630.0, -22630.0)
In [6]:
Z_new = np.copy(topo.Z)
Z_new[mask_dry_onshore] = 0.1

## plot
fig = plt.figure(figsize=(12,8))
# Z_allow_wet
ax1 = fig.add_subplot(121)
pc = pcolorcells(topo.X, topo.Y, Z_new, cmap=newcmap, vmin=crange[0], vmax=crange[1])
cbar = fig.colorbar(pc)
ax1.axis('scaled')
# diff 
ax2 = fig.add_subplot(122)
pc = pcolorcells(topo.X, topo.Y, Z_new - topo.Z, cmap=plt.get_cmap('bwr', 20), vmin=-5.0, vmax=5.0)
cbar = fig.colorbar(pc, extend='both')
ax2.axis('scaled')
Out[6]:
(-11745.0, 7785.0, -49645.0, -22615.0)
In [7]:
## write
toponew = copy.deepcopy(topo)
Z_new.fill
toponew.set_xyZ(topo.X, topo.Y, Z_new)
toponew.write(topofile_new, topo_type=3, no_data_value=topo.no_data_value, fill_value=None, header_style='geoclaw', Z_format='%15.7e', grid_registration=None)
In [8]:
## read
topo_new_check = topotools.Topography(topofile_new, topo_type=3)
In [9]:
## check 
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(121)
pc = pcolorcells(topo.X, topo.Y, topo.Z, cmap=newcmap, vmin=crange[0], vmax=crange[1])
ax.axis('scaled')
cbar = fig.colorbar(pc)
ax = fig.add_subplot(122)
pc = pcolorcells(topo_new_check.X, topo_new_check.Y, topo_new_check.Z, cmap=newcmap, vmin=crange[0], vmax=crange[1])
ax.axis('scaled')
cbar = fig.colorbar(pc)