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 = '../zone07_lonlat/depth_0030-01_zone07_lonlat.asc'
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, -50,200, 0)
fig = plt.figure()
ax = fig.add_subplot(111)
pc = pcolorcells(topo.X, topo.Y, topo.Z, cmap=newcmap, vmin=-50, vmax=200)
ax.axis('scaled')
cbar = fig.colorbar(pc)
wet_points = marching_front.select_by_flooding(topo.Z, Z1=-8.0, Z2=-1.0, 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 = pcolorcells(topo.X, topo.Y, Zdry, cmap=newcmap, vmin=-50, vmax=200)
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 = -8, Z2 = -1, max_iters=1681132 Done after 603 iterations with 109678 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,8))
# Z_allow_wet
ax1 = fig.add_subplot(121)
pc = pcolorcells(topo.X, topo.Y, Z_allow_wet, cmap=newcmap, vmin=-50, vmax=200)
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')
(136.565586679, 136.98811757699062, 34.994289679000005, 35.37330205500031)
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=-50, vmax=200)
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=plt.get_cmap('bwr', 20), vmin=-5.0, vmax=5.0)
cbar = fig.colorbar(pc, extend='both')
ax2.axis('scaled')
(136.565586679, 136.98811757699062, 34.994289679000005, 35.37330205500031)
## 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,8))
ax = fig.add_subplot(121)
pc = pcolorcells(topo.X, topo.Y, topo.Z, cmap=newcmap, vmin=-50, vmax=200)
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=-50, vmax=200)
ax.axis('scaled')
cbar = fig.colorbar(pc)