Bay-Delta SELFE Tools

grid_opt Module

Module contains routines to perform grid optimization with two options:
  1. lsqr without constraint to solve Ax=b (scipy.sparse.linalg.lsqr)
  2. minimize function 1/2*||Ax-b||^2 using the L-BFGS-B algorithm (scipy.optimize.fmin_l_bfgs_b)

x is depth at nodes with respect to a nominal reference surface (which may be greater than sea level in upstream locations).

Regularization:

  1. close to the original values for all nodes: damp
  2. minimize the 1st derivative of node elevation along land boundaries: damp_shoreline

Notes:

  • The function “cal_z_ref” only applies to the Bay Delta grid and need to be modified for other grids.
  • This script create an optimized gr3 file (_opt.gr3) only when one set of optimization parameters is specified.

First version, completed on 9/11/2013

grid_opt.cal_distance(x1, y1, x2, y2)
grid_opt.cal_length(a)
grid_opt.cal_tri_area(a)
grid_opt.cal_z_ref(xy_coor)

Define reference water surface at locations specified in xy_coor based on the assumptions, which only appliable to Bay Delta grid (1) zref increases linerly from west (ocean) to east (2) east of (x_old_river, y_old_river), zref increases lineraly towards the south

grid_opt.construct_A_and_b(damp, damp_bnd, face_coeff, vol_coeff, A_tri, A_face, A_bnd, depth_q, face_hq, bnd_h, nnode, h0_node, solver='L-BFGS-B')

construct matrix A and vector b based on specified optimization parameters damp: regularization to initial values A_bnd, damp_bnd: regularization to minimize gradient along the shoreline A_face, face_coeff: related to vertical faces of elements A_tri, vol_coeff: related to volumes of elements depth_q, face_hq: quadrature depths derived from elements and and faces bnd_h: target depth along the shoreline h0_node: original depth at nodes nnode: node number solver: either L-BFGS-B or lsqr

grid_opt.create_arg_parser()
grid_opt.create_driver()
grid_opt.create_element_shapefile(file_path, element_data, node_data, n_node_per_elm=3, n_attribute=0, attribute_def=None)

Create a shapefile representing elements as polygons Each element one is defined by a series of nodes (zero-based in element_data), where xy coordinates are provided in node_data.

Parameters :

file_path : str

Location and name of the shapefile to be created, string

element_data : array

Node number for each element

node_data : array

x,y and attribute values for each node, numpy array

n_node_per_elm : int

number of nodes (vertices) per element, deafult = 3

n_attribute : int

number of attributes defined for individual elements, default = 0

attribute_def : List

name and type of point attributes, list [[name, datatype]]

valid datatype: float (=> OFTReal), int (=> OFTInteger), str (=>OFTString) :

prototype datatype

grid_opt.create_face_shapefile(file_path, face_data, node_data, n_attribute=0, attribute_def=None)

Create a shapefile for edges

Shapefile is based on based on node # specified in face_data and xy coordinates specifed in node_data, FID is zero-based, i.e. starting at 0

Parameters :

file_path : str

location and name of the shapefile to be created, string

face_data : array

two node numbers for each face, numpy array

node_data : array

x,y and attribute values for each points, numpy array

n_attribute : int

number of attribute defined for the face feature

attribute_def : List

name and type of point attributes, list [[name, datatype]]

valid datatype: float (=> OFTReal), int (=> OFTInteger), str (=>OFTString)

grid_opt.create_node_shapefile(file_path, nodes)

Create shapefile for nodes

Parameters :

file_path : str

Location of the shapefile to be created, string

nodes : array

x,y, depth for each points, numpy array

grid_opt.create_partial_domain(elements, nodes, land_boundaries, polygon)

Select nodes, elements, and land_boundaries based on specified polygon

grid_opt.create_point_shapefile(file_path, point_data, n_attribute=0, attribute_def=None)

Create point shapefile based on xy coordinates specifed in point_data, FID is zero-based, i.e. starting at 0

Parameters :

file_path : str

location and name of the shapefile to be created, string

attribute_def : str

name and type of point attributes, list [[name, datatype]] valid datatype: float (=> OFTReal), int (=> OFTInteger), str (=>OFTString)

point_data : array

x,y and attribute values for each points, numpy array

grid_opt.create_ref_surf(elements, nodes)

create reference surface at individual nodes: step 1: derive max levation within connected elements setp 2: compare with zref calculated by cal_z_ref() and use the higher value between the two

grid_opt.depth_opt(damp, damp_bnd, face_coeff, vol_coeff, A_tri, A_face, A_bnd, depth_q, face_hq, bnd_h, nnode, h0_node, solver='L-BFGS-B')
grid_opt.fprime(x, A, b)

The gradient of the objective function with respect to the heights (x). This gradient is a vector the same size as x. A^T (Ax-b)

grid_opt.grid_opt(nodes, elements, land_boundaries, dem_list_file, opt_parms, output_files=False, outfile_prefix='', remove_boundary_faces=False, solver='L-BFGS-B')

Perform grid optimization with two solvers: (1) L-BFGS-B (default): minimize function 1/2*||Ax-b||^2 (2) linear least square solver (Ax=b). nodes: (x,y,z) for each node, numpy array, z is elevation (a place holder, valuse not used) elements: node numbers (3) for each element, numpy array land_boundaries: (land boubdary node pair index, node 1, node 2, segment no, sequence in segment), numpy array dem_list_file: a text file, containing a list of DEM used to derive elevation opt_parms: [[damp, damp_shoreline, face_coeff, volume_coeff]], list output_files: True or False, whether to create output files outfile_prefix: Prefix used for names of output files

grid_opt.main()
grid_opt.obj_func(x, A, b)

define the function to be minimize using L-BFGS-B algorithm 1/2*||Ax-b||^2

grid_opt.opt_parm_check(opt_parm_list)

check whether specified optimization parameters are valid

grid_opt.point_inside_polygon(x, y, poly)

Decide if a point is inside a polygon using the ray casting method

Parameters :

x,y : float

Point to consider

poly: List :

list of pairs (x,y) containing the coordinates of the polygon’s vertices.

Returns :

inside : boolean

True if inside

grid_opt.read_2dm(infile_2dm)
grid_opt.read_boundary_node_csv(bnd_csv_file)

Read boundary node list file Reads as (no of adjacent node pairs, 1st node no, 2nd node no, boundary segment, sequence in segment) usually generated by other scripts, can be used with .2dm or .gr3 file without boundary node information. Mainly used during the developmnet phase, kept for debugging purpose

grid_opt.read_gr3(infile_gr3)
grid_opt.read_opt_parm(opt_parm_file)
grid_opt.write_gr3(elements, nodes, outfile)

create a .gr3 file with optimized elevations, but no boundary information