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
topofile_org = '../zone06/depth_0030-01.asc'
crange = [-20.0,100.0]
topofile_new = topofile_org.replace('.asc','_masked.asc')
topofile_new = os.path.basename(topofile_new)
topo = topotools.Topography(topofile_org, topo_type=3)
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)
wet_points = marching_front.select_by_flooding(topo.Z, Z1=-6.5, Z2=-1.0, max_iters=None)
Zdry = ma.masked_array(topo.Z, wet_points)
## plot
fig = plt.figure(figsize=(12,6))
# topography with wet points masked
ax1 = fig.add_subplot(121)
pc = pcolorcells(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 = pcolorcells(topo.X, topo.Y, wet_points)
ax2.axis('scaled')
cbar = fig.colorbar(pc)
Selecting points with Z1 = -6.5, Z2 = -1, max_iters=1231200 Done after 415 iterations with 173357 points chosen
## 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,6))
# Z_allow_wet
ax1 = fig.add_subplot(121)
pc = pcolorcells(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 = pcolorcells(topo.X, topo.Y, mask_dry_onshore, vmin=0, vmax=1)
cbar = fig.colorbar(pc)
ax2.axis('scaled')
(-61715.0, -29315.0, -161115.0, -126915.0)
Z_new = np.copy(topo.Z)
Z_new[mask_dry_onshore] = 0.1
## plot
fig = plt.figure(figsize=(12,6))
# 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')
# mask_dry_onshore
ax2 = fig.add_subplot(122)
pc = pcolorcells(topo.X, topo.Y, Z_new - topo.Z, cmap='bwr', vmin=-2.0, vmax=2.0)
cbar = fig.colorbar(pc, extend='both')
ax2.axis('scaled')
(-61715.0, -29315.0, -161115.0, -126915.0)
## 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)
## read
topo_new_check = topotools.Topography(topofile_new, topo_type=3)
## check
fig = plt.figure(figsize=(12,6))
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)