.. _ruled_rectangles: Ruled Rectangles ================ **New in Version 5.7.0.** Adapted from `this notebook `_. In past versions of AMRClaw and GeoClaw, one could specify "refinement regions" as rectangles over which a `minlevel` and `maxlevel` are specified (perhaps over some time interval), as described in :ref:`refinement-regions`. Allowing only rectangles made it very easy to check each cell for inclusion in a region, but is a severe restriction -- often a number of rectangles must be used to follow a complicated coastline, for example, or else many points ended up being refined that did not require it (e.g. onshore points that never get wet in a GeoClaw simulation). We have introduced a new data structure called a "Ruled Rectangle" that is a special type of polygon for which it is also easy to check whether a given point is inside or outside, but that is much more flexible than a rectangle. It is a special case of a "ruled surface", which is a surface in 3D that is bounded by two curves, each parameterized by :math:`s` from :math:`s_1` to :math:`s_2`, and with the surface defined as the union of all the line segments connecting one curve to the other for each value of :math:`s` in :math:`[s_1,s_2]`. A Ruled Rectangle is a special case in which each curve lies in the :math:`x`-:math:`y` plane and either :math:`s=x` for some range of :math:`x` values or :math:`s=y` for some range of :math:`y` values. If :math:`s=x`, for example, then the line segments defining the surface are intervals :math:`y_{\scriptstyle lower}(x) \leq y \leq y_{\scriptstyle upper}(x)` for each :math:`x` over some range :math:`x_1 \leq x \leq x_2`. It is easy to check if a given :math:`(x_c,y_c)` is in this region: it is if :math:`x_1 \leq x_c \leq x_2` and in addition :math:`y_{\scriptstyle lower}(x_c) \leq y_c \leq y_{\scriptstyle upper}(x_c)`. The class ``clawpack.amrclaw.region_tools.RuledRectangle`` supports a subset of ruled rectangles defined by a finite set of :math:`s` values along a coordinate line, e.g. :math:`s[0] < s[1] < s[2] < \cdots < s[N]` and for each :math:`s[k]` two values `lower[k]` and `upper[k]`. If `rr` is an instance of this class then `rr.s`, `rr.lower`, and `rr.upper` contain these arrays. Whether `s` corresponds to `x` or `y` is determined by: - If `rr.ixy in [1, 'x']` then `s` gives a set of :math:`x` values, - If `rr.ixy in [2, 'y']` then `s` gives a set of :math:`y` values. The points specified can then be connected by line segments to define a Ruled Rectangle, and this is done if `rr.method == 1` (piecewise linear). On the other hand, if `rr.method == 0` (piecewise constant) then the values `lower[k], upper[k]` are used as the bounds for all `s` in the interval :math:`s[k] \leq s \leq s[k+1]` (for :math:`k = 0,~1,~\ldots,~N-1`). In this case the values `lower[N], upper[N]` are not used. This also defines a polygon, but one that consists of a set of *stacked boxes*. The advantage of the latter form is that it is slightly easier to check if a point is in the Ruled Rectangle since no linear interpolation is required along the edges. Also for some applications we want the Ruled Rectangle to exactly cover a contiguous set of finite volume grid cells, which has the shape of a set of stacked boxes. Contents -------- - `Examples <#rr-examples>`__ - `Relation to convex polygons <#rr-polygons>`__ - `Other attributes and methods <#rr-attr>`__ - `A simple rectangle <#rr-rect>`__ - `Defining a Ruled Rectangle covering selected cells <#rr-cover>`__ .. _rr-examples: Examples -------- Some simple examples follow: Define a Ruled Rectangle by specifying a set of points: .. code:: ipython3 from clawpack.amrclaw import region_tools rr = region_tools.RuledRectangle() rr.ixy = 1 # so s refers to x, lower & upper are limits in y rr.s = array([1,2,4,6]) rr.lower = array([0,2,1,3]) rr.upper = array([4,5,3,6]) Setting `rr.method` to 1 or 0 gives a Ruled Rectangle in which the points specified above are connected by lines or used to define stacked boxes. Both are illustrated in the figure below. Note that we use the method `rr.vertices()` to return a list of all the vertices of the polygon defined by `rr` for plotting purposes. .. code:: ipython3 figure(figsize=(10,5)) subplot(121) rr.method = 1 xv,yv = rr.vertices() plot(xv,yv,'b') grid(True) axis([0,7,-1,6]) title('With method==1') subplot(122) rr.method = 0 xv,yv = rr.vertices() plot(xv,yv,'r') grid(True) axis([0,7,-1,6]) title('With method==0') .. image:: RuledRectangles_files/RuledRectangles_10_1.png In the plots above the `s` values correspond to `x = 1, 2, 4, 6`, and the `lower` and `upper` arrays define ranges in `y`. If we set `rr.ixy = 2` or `'y'`, then the `s` values will instead correspond to `y = 1, 2, 4, 6` and the `lower` and `upper` will define ranges in `x`. This is illustrated in the plots below. .. code:: ipython3 rr.ixy = 2 # so s refers to y, lower & upper are limits in x figure(figsize=(10,5)) subplot(121) rr.method = 1 xv,yv = rr.vertices() plot(xv,yv,'b') grid(True) axis([-1,6,0,7]) title('With method==1') subplot(122) rr.method = 0 xv,yv = rr.vertices() plot(xv,yv,'r') grid(True) axis([-1,6,0,7]) title('With method==0') .. image:: RuledRectangles_files/RuledRectangles_12_1.png .. _rr-polygons: Relation to convex polygons --------------------------- Note that the polygons above are not convex, but clearly some Ruled Rectangles would be convex. Conversely, *any* convex polygon can be expressed as a Ruled Rectangle --- simply order the vertices so that the :math:`x` values are increasing, for example, and use these as the `s` values in a `RuledRectangle` with `ixy='x'`. Then for each :math:`x` there is a connected interval of :math:`y` values that lie within the polygon (by convexity), so this defines the `lower` and `upper` values. (Or you could start by ordering vertices by increasing :math:`y` values and similarly define a `RuledRectangle` with `ixy='y'`.) So a `RuledRectangle` is a nice generalization of a convex polygon for which it is easy to check inclusion of an arbitrary point. .. _rr-attr: Other attributes and methods ---------------------------- `ds` ~~~~~~ If the points `s[k]` are equally spaced then `ds` is the spacing between them. This makes it quicker to determine what two points an arbitrary value of :math:`s` lies between when determining whether a large set of points are inside or outside the Ruled Rectangle, rather than having to search. The Ruled Rectangle defined above has unequally spaced points and the `ds` attribute is set to `-1` in this case. .. code:: ipython3 rr.ds .. parsed-literal:: -1 `slu` ~~~~~~~ Rather than specifying `s`, `lower`, and `upper` separately, you can specify an array `slu` with three columns in defining a `RuledRectangle`, and such an array is returned by the `slu` method: .. code:: ipython3 rr.slu() .. parsed-literal:: array([[1, 0, 4], [2, 2, 5], [4, 1, 3], [6, 3, 6]]) Here's an example defining a `RuledRectangle` via `slu`: .. code:: ipython3 slu = vstack((linspace(0,14,8), zeros(8), [1,2,1,2,1,2,1,2])).T print('slu = \n', slu) rr = region_tools.RuledRectangle(slu=slu) rr.ixy = 2 rr.method = 1 xv,yv = rr.vertices() plot(xv,yv) .. parsed-literal:: slu = [[ 0. 0. 1.] [ 2. 0. 2.] [ 4. 0. 1.] [ 6. 0. 2.] [ 8. 0. 1.] [10. 0. 2.] [12. 0. 1.] [14. 0. 2.]] .. image:: RuledRectangles_files/RuledRectangles_19_2.png `bounding_box` ~~~~~~~~~~~~~~~~ `rr.bounding_box()` returns the smallest rectangle `[x1,x2,y1,y2]` containing the ruled rectangle: .. code:: ipython3 rr.bounding_box() .. parsed-literal:: [0.0, 2.0, 0.0, 14.0] `mask_outside` ~~~~~~~~~~~~~~~~ If `X,Y` are 2D numpy arrays defining `(x,y)` coordinates on a grid, then `rr.mask_outside(X,Y)` returns a mask array `M` of the same shape as `X,Y` that is `True` at points outside the Ruled Rectangle. .. code:: ipython3 x = linspace(0,3,31) y = linspace(0,16,81) X,Y = meshgrid(x,y) Z = X + Y # sample data values to plot M = rr.mask_outside(X,Y) Zm = ma.masked_array(Z,mask=M) pcolorcells(X,Y,Zm) colorbar() .. image:: RuledRectangles_files/RuledRectangles_23_1.png read and write, and instantiating from a file ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ `rr.write(fname)` writes out the `slu` array and other attributes to file `fname`, and `rr.read(fname)` can be used to read in such a file. You can also specify `fname` when instantiating a new Ruled Rectangle: .. code:: ipython3 rr.write('RRzigzag.data') .. code:: ipython3 rr2 = region_tools.RuledRectangle('RRzigzag.data') rr2.slu() .. parsed-literal:: array([[ 0., 0., 1.], [ 2., 0., 2.], [ 4., 0., 1.], [ 6., 0., 2.], [ 8., 0., 1.], [10., 0., 2.], [12., 0., 1.], [14., 0., 2.]]) Here's what the file looks like: .. code:: ipython3 lines = open('RRzigzag.data').readlines() for line in lines: print(line.strip()) .. parsed-literal:: 2 ixy 1 method 2 ds 8 nrules 0.000000000 0.000000000 1.000000000 2.000000000 0.000000000 2.000000000 4.000000000 0.000000000 1.000000000 6.000000000 0.000000000 2.000000000 8.000000000 0.000000000 1.000000000 10.000000000 0.000000000 2.000000000 12.000000000 0.000000000 1.000000000 14.000000000 0.000000000 2.000000000 Note that this Ruled Rectangle has equally spaced points and so `ds = 2` is the spacing. `make_kml` ~~~~~~~~~~~~ `rr.make_kml()` can be used to create a kml file that can be opened on Google Earth to show the polygon defined by `rr`. This assumes that `x` corresponds to longitude and `y` to latitude and is designed for GeoClaw applications. Several optional arguments can be specified: `fname, name, color, width, verbose`. A GeoClaw AMR flag region ~~~~~~~~~~~~~~~~~~~~~~~~~ The figure below shows a Ruled Rectangle designed to cover Admiralty Inlet, the water way between the Kitsap Peninsula and Whidbey Island connecting the Strait of Juan de Fuca to lower Puget Sound. For some tsunami modeling problems it is important to cover this region with a finer grid than is needed elsewhere. .. code:: ipython3 Image('figs/RuledRectangle_AdmiraltyInlet.png', width=400) .. image:: RuledRectangles_files/RuledRectangles_32_0.png :width: 400px The Ruled Rectangle shown above was defined by the code below: .. code:: ipython3 slu = \ array([[ 47.851, -122.75 , -122.300], [ 47.955, -122.75 , -122.300], [ 48. , -122.8 , -122.529], [ 48.036, -122.8 , -122.578], [ 48.12 , -122.9 , -122.577], [ 48.187, -122.9 , -122.623], [ 48.191, -122.9 , -122.684], [ 48.221, -122.9 , -122.755]]) rr_admiralty = region_tools.RuledRectangle(slu=slu) rr_admiralty.ixy = 'y' rr_admiralty.method = 1 rr_name = 'RuledRectangle_AdmiraltyInlet' rr_admiralty.write(rr_name + '.data') rr_admiralty.make_kml(fname=rr_name+'.kml', name=rr_name) The file `RuledRectangle_AdmiraltyInlet.data` can then be used as a "flag region" in the modified GeoClaw code, see `FlagRegions.ipynb `__ for more details. The file `RuledRectangle_AdmiraltyInlet.kml` can be opened in Google Earth to show the polygon, as captured in the figure above. .. _rr-rect: A simple rectangle ------------------ A simple rectangle with extent `[x1,x2,y1,y2]` can be specified as a `RuledRectangle` via e.g. : :: rectangle = region_tools.RuledRectangle() rectangle.ixy = 'x' rectangle.s = [x1, x2] rectangle.lower = [y1, y1] rectangle.upper = [y2, y2] rectangle.method = 0 This can be done for you when instantiating a `RuledRectangle` using: :: rectangle = region_tools.RuledRectangle(rect=[x1,x2,y1,y2]) For example: .. code:: ipython3 rectangle = region_tools.RuledRectangle(rect=[1,3,5,6]) xv,yv = rectangle.vertices() plot(xv,yv,'b') grid(True) axis([0,4,4,7]); .. image:: RuledRectangles_files/RuledRectangles_37_0.png .. _rr-cover: Defining a Ruled Rectangle covering selected cells -------------------------------------------------- The module function ``region_tools.ruledrectangle_covering_selected_points`` can be used to generate a Ruled Rectangle that covers a specified set of points as compactly as possible. This is useful for generating AMR refinement regions that cover a set of points where we want to enforce a fine grid without including too many other points. First generate some sample data: .. code:: ipython3 x_edge = linspace(-5,27,33) y_edge = linspace(-7,7,15) x_center = 0.5*(x_edge[:-1] + x_edge[1:]) y_center = 0.5*(y_edge[:-1] + y_edge[1:]) X,Y = meshgrid(x_center,y_center) pts_chosen = where(abs(X-Y**2) < 4., 1, 0) pts_chosen = where(logical_or(X>24, Y<-4), 0, pts_chosen) pcolorcells(x_edge,y_edge,pts_chosen, cmap=cm.Blues, edgecolor=[.8,.8,.8], linewidth=0.1) .. image:: RuledRectangles_files/RuledRectangles_39_1.png The array `pts_chosen` has been defined with the value 1 in the dark blue cells in the figure above, and 0 elsewhere. In this case the points can be covered by a Ruled Rectangle with `ixy = 'y'` most efficiently, giving a polygon that covers only the selected points. Note we use `method = 0` to generate a Ruled Rectangle that covers all the grid cells: .. code:: ipython3 rr = region_tools.ruledrectangle_covering_selected_points(x_center, y_center, pts_chosen, ixy='y', method=0, verbose=True) pcolorcells(x_edge,y_edge,pts_chosen, cmap=cm.Blues, edgecolor=[.8,.8,.8], linewidth=0.1) xv,yv = rr.vertices() plot(xv,yv,'r',label='Ruled Rectangle') legend(loc='lower right') title('With ixy=2 and method=0'); .. parsed-literal:: Extending rectangles to cover grid cells RuledRectangle covers 80 grid points .. image:: RuledRectangles_files/RuledRectangles_42_1.png By contrast, if we use `ixy = 'x'`, the minimal Ruled Rectangle covering the selected cells will also cover a number of cells that were not selected: .. code:: ipython3 rr = region_tools.ruledrectangle_covering_selected_points(x_center, y_center, pts_chosen, ixy='x', method=0, verbose=True) pcolorcells(x_edge,y_edge,pts_chosen, cmap=cm.Blues, edgecolor=[.8,.8,.8], linewidth=0.1) xv,yv = rr.vertices() plot(xv,yv,'r',label='Ruled Rectangle') legend(loc='lower right') title('With ixy=1 and method=0'); .. parsed-literal:: Extending rectangles to cover grid cells RuledRectangle covers 129 grid points .. image:: RuledRectangles_files/RuledRectangles_44_1.png Note that if we use `method = 1` then the Ruled Rectangle covers the center of each cell but not the entire grid cell for cells near the boundary: .. code:: ipython3 rr = region_tools.ruledrectangle_covering_selected_points(x_center, y_center, pts_chosen, ixy='y', method=1, verbose=True) pcolorcells(x_edge,y_edge,pts_chosen, cmap=cm.Blues, edgecolor=[.8,.8,.8], linewidth=0.1) xv,yv = rr.vertices() plot(xv,yv,'r',label='Ruled Rectangle') legend(loc='lower right') title('With ixy=2 and method=1'); .. parsed-literal:: RuledRectangle covers 63 grid points .. image:: RuledRectangles_files/RuledRectangles_46_1.png With `ixy='x'` and `method=1` the Ruled Rectangle degenerates in the upper right corner to a line segment that covers only the cell centers: .. code:: ipython3 rr = region_tools.ruledrectangle_covering_selected_points(x_center, y_center, pts_chosen, ixy='x', method=1, verbose=True) pcolorcells(x_edge,y_edge,pts_chosen, cmap=cm.Blues, edgecolor=[.8,.8,.8], linewidth=0.1) xv,yv = rr.vertices() plot(xv,yv,'r',label='Ruled Rectangle') legend(loc='lower right') title('With ixy=1 and method=1'); .. parsed-literal:: RuledRectangle covers 100 grid points .. image:: RuledRectangles_files/RuledRectangles_48_1.png Example covering the continental shelf ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The figure below shows a Ruled Rectangle that roughly covers the continental shelf offshore of Vancouver Island and Washington state. The region may need to be refined to a higher level than the deeper ocean in order to capture shoaling tsunami waves interacting with the shelf topography. This region was defined by first selecting a set of points from etopo1 topography satisfying certain constraints on elevation using the marching front algorithm described in :ref:`marching_front`, and then using the ``region_tools.ruledrectangle_covering_selected_points`` function to build a Ruled Rectangle covering these points. .. image:: RuledRectangles_files/RuledRectangles_50_0.png :width: 400px