agg_wy2013/dwaq_aggregate_debug.py at master · rustychris/agg_wy2013 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
"""
Sample driver script for aggregation of dwaq hydrodynamics -
  this version is for debugging disappearing BCs.
"""

import matplotlib
import os
import datetime
import stompy.model.delft.waq_scenario as waq
import six

# Inputs:
agg_grid_shp="Agg2Shapefile_edited/agg_shapefile_edited.shp"
hyd_fn="/hpcvol1/public/sfb_dfm_v2/wy2013c/DFM_DELWAQ_wy2013c_adj/wy2013c.hyd"
output_fn="wy2013c-agg-dbg/wy2013c-agg.hyd"


#----------

# Clean inputs
def clean_shapefile(shp_in):
    """
    break multipolygons into individual polygons.
    """
    from stompy.spatial import wkb2shp
    geoms=wkb2shp.shp2geom(agg_grid_shp)

    multi_count=0

    new_geoms=[]
    for fi,feat in enumerate(geoms):
        if feat['geom'].type=='Polygon':
            new_geoms.append(feat['geom'])
        else:
            multi_count+=1
            for g in feat['geom'].geoms:
                new_geoms.append(g)
    if multi_count:
        cleaned=agg_grid_shp.replace('.shp','-cleaned.shp')
        assert cleaned!=agg_grid_shp
        wkb2shp.wkb2shp(cleaned,new_geoms,overwrite=True)

        return cleaned
    else:
        return shp_in


# Processing:

# remove multipolygons from inputs
shp=clean_shapefile(agg_grid_shp)

##
six.moves.reload_module(waq)

# open the original dwaq hydrodynamics
hydro_orig=waq.HydroFiles(hyd_fn)

# create object representing aggregated hydrodynamics
# sparse_layers: for z-layer inputs this can be True, in which cases cells are only output for the
#    layers in which they are above the bed.  Usually a bad idea.  Parts of DWAQ assume
#    each 2D cell exists across all layers
# agg_boundaries: if True, multiple boundary inputs entering a single aggregated cell will be
#   merged into a single boundary input.  Generally best to keep this as False.
hydro_agg=waq.HydroAggregator(hydro_in=hydro_orig,
                              agg_shp=shp,
                              sparse_layers=False,
                              agg_boundaries=False)

##

bdefs=hydro_agg.boundary_defs()

# this does *not* have burlingame
# probably related to how burlingame has the exact same x values as millbrae, i.e.
# it goes into the same segment in the original grid.
# at first burlingame was suppressed even in hydro_orig.
# that has been fixed (by distinguishing between boundary exchanges that enter
# the same segment, and assuming a fixed order of them.)

##

gbl=hydro_agg.group_boundary_links()

gbl_unagg=hydro_orig.group_boundary_links()
##

# original code finds 4 links here, 3 real, 1 boundary.
# new code finds the extra boundary link here

##


# The code to write dwaq hydro is wrapped up in the code to write a dwaq model inp file,
# so we pretend to set up a dwaq simulation, even though the goal is just to write
# the hydro.
name=os.path.basename(output_fn.replace('.hyd',''))
class Writer(waq.Scenario):
    name=name
    desc=(name,
          agg_grid_shp,
          'aggregated')
    # output directory inferred from output hyd path
    base_path=os.path.dirname(output_fn)

# Define the subset of timesteps to write out, in this case the
# whole run.
sec=datetime.timedelta(seconds=1)
start_time=hydro_agg.time0+hydro_agg.t_secs[ 0]*sec
stop_time =hydro_agg.time0+hydro_agg.t_secs[-1]*sec

# probably would have been better to just pass name, desc, base_path in here,
# rather than using a shell subclass.
writer=Writer(hydro=hydro_agg,
              start_time=start_time,
              stop_time=stop_time)

# This step is super slow.  Watch the output directory for progress.
# Takes ~20 hours on HPC for the full wy2013 run.
writer.cmd_write_hydro()