Identify storms in gridded model output and extract uniform sized patches around the storm centers of mass.

        Returns:
196703
def find_model_patch_tracks(self):
        
        self.model_grid.load_data()
        tracked_model_objects = []
        model_objects = []
        if self.model_grid.data is None:
            print("No model output found")
            return tracked_model_objects
        min_orig = self.model_ew.min_thresh
        max_orig = self.model_ew.max_thresh
        data_increment_orig = self.model_ew.data_increment
        self.model_ew.min_thresh = 0
        self.model_ew.data_increment = 1
        self.model_ew.max_thresh = 100
        for h, hour in enumerate(self.hours):
            # Identify storms at each time step and apply size filter
            print("Finding {0} objects for run {1} Hour: {2:02d}".format(self.ensemble_member,
                                                                         self.run_date.strftime("%Y%m%d%H"), hour))
            if self.mask is not None:
                model_data = self.model_grid.data[h] * self.mask
            else:
                model_data = self.model_grid.data[h]
            model_data[:self.patch_radius] = 0
            model_data[-self.patch_radius:] = 0
            model_data[:, :self.patch_radius] = 0
            model_data[:, -self.patch_radius:] = 0
            scaled_data = np.array(rescale_data(model_data, min_orig, max_orig))
            hour_labels = label_storm_objects(scaled_data, "ew",
                                              self.model_ew.min_thresh, self.model_ew.max_thresh,
                                              min_area=self.size_filter, max_area=self.model_ew.max_size,
                                              max_range=self.model_ew.delta, increment=self.model_ew.data_increment,
                                              gaussian_sd=self.gaussian_window)
            model_objects.extend(extract_storm_patches(hour_labels, model_data, self.model_grid.x,
                                                       self.model_grid.y, [hour],
                                                       dx=self.model_grid.dx,
                                                       patch_radius=self.patch_radius))
            for model_obj in model_objects[-1]:
                dims = model_obj.timesteps[-1].shape
                if h > 0:
                    model_obj.estimate_motion(hour, self.model_grid.data[h-1], dims[1], dims[0])
            del scaled_data
            del model_data
            del hour_labels
        tracked_model_objects.extend(track_storms(model_objects, self.hours,
                                                  self.object_matcher.cost_function_components,
                                                  self.object_matcher.max_values,
                                                  self.object_matcher.weights))
        self.model_ew.min_thresh = min_orig
        self.model_ew.max_thresh = max_orig
        self.model_ew.data_increment = data_increment_orig
        return tracked_model_objects
130708
def get_xy_grids(ds, stride=1, getval=False):
    
    gt = ds.GetGeoTransform()
    #stride = stride_m/gt[1]
    pX = np.arange(0, ds.RasterXSize, stride)
    pY = np.arange(0, ds.RasterYSize, stride)
    psamp = np.meshgrid(pX, pY)
    mX, mY = pixelToMap(psamp[0], psamp[1], gt)
    return mX, mY
1154865
def _compute_grid(self):
        
        widths = []
        for c in range(self.cols):
            c_widths = []
            for r in range(self.rows):
                subplot = self.subplots.get((r, c), None)
                nsubplots = 1 if subplot is None else len(subplot.layout)
                c_widths.append(2 if nsubplots > 1 else 1)
            widths.append(max(c_widths))

        heights = []
        for r in range(self.rows):
            r_heights = []
            for c in range(self.cols):
                subplot = self.subplots.get((r, c), None)
                nsubplots = 1 if subplot is None else len(subplot.layout)
                r_heights.append(2 if nsubplots > 2 else 1)
            heights.append(max(r_heights))

        # Generate empty grid
        rows = sum(heights)
        cols = sum(widths)
        grid = [[None]*cols for _ in range(rows)]

        return grid
730628
def grid_list(data):
    '''
    
    '''
    a = []
    b = []
    M = []
    for i in data:
        a.append(i[0])  # np.array([i[1] for i in data], dtype=np.float64)
        b.append(i[1])  # np.array([i[0] for i in data], dtype=np.float64)
        M.append(i[2])  # np.array([i[2] for i in data], dtype=np.float64)
    a = np.array(a, dtype=np.float64).tolist()
    b = np.array(b, dtype=np.float64).tolist()
    M = np.array(M, dtype=np.float64).tolist()
    a = list(set(a))
    b = list(set(b))
    return a, b, M
559254
def rain_series(self):
        
        return [(tstamp, \
                self._station_history.get_measurements()[tstamp]['rain']) \
                for tstamp in self._station_history.get_measurements()]
734600
def predict_size_models(self, model_names,
                            input_columns,
                            metadata_cols,
                            data_mode="forecast"):
        
        groups = self.size_models.keys()
        predictions = {}
        for group in groups:
            group_data = self.data[data_mode]["combo"].loc[self.data[data_mode]["combo"][self.group_col] == group]
            if group_data.shape[0] > 0:
                predictions[group] = {}
                output_values = self.size_models[group]["outputvalues"].astype(int)
                for m, model_name in enumerate(model_names):
                    print("{0} {1}".format(group, model_name))
                    pred_col_names = [model_name.replace(" ", "-") + "_{0:02d}".format(p) for p in output_values]
                    predictions[group][model_name] = group_data[metadata_cols]
                    pred_vals = self.size_models[group][model_name].predict_proba(group_data[input_columns])
                    pred_classes = self.size_models[group][model_name].classes_
                    pred_pdf = np.zeros((pred_vals.shape[0], output_values.size))
                    for pcv, pc in enumerate(pred_classes):
                        idx = np.where(output_values == pc)[0][0]
                        pred_pdf[:, idx] = pred_vals[:, pcv]
                    for pcn, pred_col_name in enumerate(pred_col_names):
                        predictions[group][model_name].loc[:, pred_col_name] = pred_pdf[:, pcn]
        return predictions
695609
def gridPointsFromEdges(edges, nCells, dtype=float):
    
    e = edges
    sx, sy = nCells[0] + 1, nCells[1] + 1
    # horizontal lines
    x0 = np.linspace(e[0, 0], e[3, 0], sy, dtype=dtype)
    x1 = np.linspace(e[1, 0], e[2, 0], sy, dtype=dtype)

    y0 = np.linspace(e[0, 1], e[3, 1], sy, dtype=dtype)
    y1 = np.linspace(e[1, 1], e[2, 1], sy, dtype=dtype)
    # points:
    p = np.empty(shape=(sx * sy, 2))
    n0 = 0
    n1 = sx
    for x0i, x1i, y0i, y1i in zip(x0, x1, y0, y1):
        p[n0:n1, 0] = np.linspace(x0i, x1i, sx)
        p[n0:n1, 1] = np.linspace(y0i, y1i, sx)
        n0 = n1
        n1 += sx
    return p
782404
def get_surface_boundaries(self):
        
        min_lon, min_lat, max_lon, max_lat = self.get_bounding_box()
        return [[min_lon, max_lon]], [[min_lat, max_lat]]
979496
def _parse_coords(coord_lines):
        
        paras = {}
        var_pattern = re.compile(r"^([A-Za-z]+\S*)[\s=,]+([\d\-\.]+)$")
        for l in coord_lines:
            m = var_pattern.match(l.strip())
            if m:
                paras[m.group(1).strip("=")] = float(m.group(2))

        species = []
        coords = []
        # Stores whether a Zmatrix format is detected. Once a zmatrix format
        # is detected, it is assumed for the remaining of the parsing.
        zmode = False
        for l in coord_lines:
            l = l.strip()
            if not l:
                break
            if (not zmode) and GaussianInput._xyz_patt.match(l):
                m = GaussianInput._xyz_patt.match(l)
                species.append(m.group(1))
                toks = re.split(r"[,\s]+", l.strip())
                if len(toks) > 4:
                    coords.append([float(i) for i in toks[2:5]])
                else:
                    coords.append([float(i) for i in toks[1:4]])
            elif GaussianInput._zmat_patt.match(l):
                zmode = True
                toks = re.split(r"[,\s]+", l.strip())
                species.append(toks[0])
                toks.pop(0)
                if len(toks) == 0:
                    coords.append(np.array([0, 0, 0]))
                else:
                    nn = []
                    parameters = []
                    while len(toks) > 1:
                        ind = toks.pop(0)
                        data = toks.pop(0)
                        try:
                            nn.append(int(ind))
                        except ValueError:
                            nn.append(species.index(ind) + 1)
                        try:
                            val = float(data)
                            parameters.append(val)
                        except ValueError:
                            if data.startswith("-"):
                                parameters.append(-paras[data[1:]])
                            else:
                                parameters.append(paras[data])
                    if len(nn) == 1:
                        coords.append(np.array([0, 0, parameters[0]]))
                    elif len(nn) == 2:
                        coords1 = coords[nn[0] - 1]
                        coords2 = coords[nn[1] - 1]
                        bl = parameters[0]
                        angle = parameters[1]
                        axis = [0, 1, 0]
                        op = SymmOp.from_origin_axis_angle(coords1, axis,
                                                           angle, False)
                        coord = op.operate(coords2)
                        vec = coord - coords1
                        coord = vec * bl / np.linalg.norm(vec) + coords1
                        coords.append(coord)
                    elif len(nn) == 3:
                        coords1 = coords[nn[0] - 1]
                        coords2 = coords[nn[1] - 1]
                        coords3 = coords[nn[2] - 1]
                        bl = parameters[0]
                        angle = parameters[1]
                        dih = parameters[2]
                        v1 = coords3 - coords2
                        v2 = coords1 - coords2
                        axis = np.cross(v1, v2)
                        op = SymmOp.from_origin_axis_angle(
                            coords1, axis, angle, False)
                        coord = op.operate(coords2)
                        v1 = coord - coords1
                        v2 = coords1 - coords2
                        v3 = np.cross(v1, v2)
                        adj = get_angle(v3, axis)
                        axis = coords1 - coords2
                        op = SymmOp.from_origin_axis_angle(
                            coords1, axis, dih - adj, False)
                        coord = op.operate(coord)
                        vec = coord - coords1
                        coord = vec * bl / np.linalg.norm(vec) + coords1
                        coords.append(coord)

        def _parse_species(sp_str):
            """
            The species specification can take many forms. E.g.,
            simple integers representing atomic numbers ("8"),
            actual species string ("C") or a labelled species ("C1").
            Sometimes, the species string is also not properly capitalized,
            e.g, ("c1"). This method should take care of these known formats.
            """
            try:
                return int(sp_str)
            except ValueError:
                sp = re.sub(r"\d", "", sp_str)
                return sp.capitalize()

        species = [_parse_species(sp) for sp in species]

        return Molecule(species, coords)
1123089
def wind_series(self):
        
        return [(timestamp, \
                self._station_history.get_measurements()[timestamp]['wind']) \
                for timestamp in self._station_history.get_measurements()]
112541
def get_window_forecasts(self):
        
        for model_name in self.model_names:
            self.window_forecasts[model_name] = {}
            for size_threshold in self.size_thresholds:
                self.window_forecasts[model_name][size_threshold] = \
                    np.array([self.raw_forecasts[model_name][size_threshold][sl].sum(axis=0)
                              for sl in self.hour_windows])
530892
def _get_flat_core_sizes(cores):
  
  core_sizes_lists = []
  for core in cores:
    flat_output_size = nest.flatten(core.output_size)
    core_sizes_lists.append(
        [tf.TensorShape(size).as_list() for size in flat_output_size])
  return core_sizes_lists
386732
def compute_eigh(self):
        
        self.eigh1 = _SHGrid.from_array(_np.zeros_like(self.vxx.data),
                                        grid='DH')
        self.eigh2 = _SHGrid.from_array(_np.zeros_like(self.vxx.data),
                                        grid='DH')
        self.eighh = _SHGrid.from_array(_np.zeros_like(self.vxx.data),
                                        grid='DH')

        for i in range(self.nlat):
            for j in range(self.nlon):
                a = _np.array([[self.vxx.data[i, j],
                                self.vxy.data[i, j]],
                               [self.vyx.data[i, j],
                                self.vyy.data[i, j]]])

                eigs = _eigvalsh(a)

                self.eigh1.data[i, j] = eigs[1]
                self.eigh2.data[i, j] = eigs[0]

                if abs(eigs[0]) >= abs(eigs[1]):
                    self.eighh.data[i, j] = eigs[0]
                else:
                    self.eighh.data[i, j] = eigs[1]
971385
def _initialize_model_diagnostics(self):
        

        self.add_diagnostic('entspec',
                description='barotropic enstrophy spectrum',
                function= (lambda self:
                    np.abs((self.Hi[:,np.newaxis,np.newaxis]*self.qh).sum(axis=0))**2/self.H) )

        self.add_diagnostic('KEspec_modal',
                description='modal KE spectra',
                function= (lambda self:
                    self.wv2*(np.abs(self.phn)**2)/self.M**2 ))

        self.add_diagnostic('PEspec_modal',
                description='modal PE spectra',
                function= (lambda self:
                    self.kdi2[1:,np.newaxis,np.newaxis]*(np.abs(self.phn[1:,:,:])**2)/self.M**2 ))

        self.add_diagnostic('APEspec',
                description='available potential energy spectrum',
                function= (lambda self:
                           (self.f2gpi*
                            np.abs(self.ph[:-1]-self.ph[1:])**2).sum(axis=0)/self.H))

        self.add_diagnostic('KEflux',
                    description='spectral divergence of flux of kinetic energy',
                    function =(lambda self: (self.Hi[:,np.newaxis,np.newaxis]*
                               (self.ph.conj()*self.Jpxi).real).sum(axis=0)/self.H))

        self.add_diagnostic('APEflux',
                    description='spectral divergence of flux of available potential energy',
                    function =(lambda self: (self.Hi[:,np.newaxis,np.newaxis]*
                               (self.ph.conj()*self.JSp).real).sum(axis=0)/self.H))

        self.add_diagnostic('APEgenspec',
                    description='the spectrum of the rate of generation of available potential energy',
                    function =(lambda self: (self.Hi[:,np.newaxis,np.newaxis]*
                                (self.Ubg[:,np.newaxis,np.newaxis]*self.k +
                                 self.Vbg[:,np.newaxis,np.newaxis]*self.l)*
                                (1j*self.ph.conj()*self.Sph).real).sum(axis=0)/self.H))

        self.add_diagnostic('ENSflux',
                 description='barotropic enstrophy flux',
                 function = (lambda self: (-self.Hi[:,np.newaxis,np.newaxis]*
                              (self.qh.conj()*self.Jq).real).sum(axis=0)/self.H))

        self.add_diagnostic('ENSgenspec',
                    description='the spectrum of the rate of generation of barotropic enstrophy',
                    function = (lambda self:
                            -(self.Hi[:,np.newaxis,np.newaxis]*((self.ikQy -
                            self.ilQx)*(self.Sph.conj()*self.ph)).real).sum(axis=0)/self.H))
879739
def iterGrid(self, minZoom, maxZoom):
        ""
        assert minZoom in range(0, len(self.RESOLUTIONS))
        assert maxZoom in range(0, len(self.RESOLUTIONS))
        assert minZoom <= maxZoom

        for zoom in xrange(minZoom, maxZoom + 1):
            [minRow, minCol, maxRow, maxCol] = self.getExtentAddress(zoom)
            for row in xrange(minRow, maxRow + 1):
                for col in xrange(minCol, maxCol + 1):
                    tileBounds = self.tileBounds(zoom, col, row)
                    yield (tileBounds, zoom, col, row)
601037
def items(cls):
        
        return [
            cls.PRECIPITATION,
            cls.WIND,
            cls.TEMPERATURE,
            cls.PRESSURE
        ]
68363
def get_surface_boundaries(self):
        
        mesh = self.mesh
        lons = numpy.concatenate((mesh.lons[0, :],
                                  mesh.lons[1:, -1],
                                  mesh.lons[-1, :-1][::-1],
                                  mesh.lons[:-1, 0][::-1]))
        lats = numpy.concatenate((mesh.lats[0, :],
                                  mesh.lats[1:, -1],
                                  mesh.lats[-1, :-1][::-1],
                                  mesh.lats[:-1, 0][::-1]))
        return [lons], [lats]