calculate kriging factors (weights) for a structured grid.

        Parameters
        ----------
        spatial_reference : (flopy.utils.reference.SpatialReference)
            a spatial reference that describes the orientation and
            spatail projection of the the structured grid
        zone_array : (numpy.ndarray)
            an integer array of zones to use for kriging.  If not None,
            then point_data must also contain a "zone" column.  point_data
            entries with a zone value not found in zone_array will be skipped.
            If None, then all point_data will (potentially) be used for
            interpolating each grid node. Default is None
        minpts_interp : (int)
            minimum number of point_data entires to use for interpolation at
            a given grid node.  grid nodes with less than minpts_interp
            point_data found will be skipped (assigned np.NaN).  Defaut is 1
        maxpts_interp : (int)
            maximum number of point_data entries to use for interpolation at
            a given grid node.  A larger maxpts_interp will yield "smoother"
            interplation, but using a large maxpts_interp will slow the
            (already) slow kriging solution process and may lead to
            memory errors. Default is 20.
        search_radius : (float)
            the size of the region around a given grid node to search for
            point_data entries. Default is 1.0e+10
        verbose : (boolean)
            a flag to  echo process to stdout during the interpolatino process.
            Default is False
        var_filename : (str)
            a filename to save the kriging variance for each interpolated grid node.
            Default is None.
        forgive : (boolean)
            flag to continue if inversion of the kriging matrix failes at one or more
            grid nodes.  Inversion usually fails if the kriging matrix is singular,
            resulting from point_data entries closer than EPSILON distance.  If True,
            warnings are issued for each failed inversion.  If False, an exception
            is raised for failed matrix inversion.

        Returns
        -------
        df : pandas.DataFrame
            a dataframe with information summarizing the ordinary kriging
            process for each grid node

        Note
        ----
        this method calls OrdinaryKrige.calc_factors()


        Example
        -------
        ``>>>import flopy``

        ``>>>import pyemu``

        ``>>>v = pyemu.utils.geostats.ExpVario(a=1000,contribution=1.0)``

        ``>>>gs = pyemu.utils.geostats.GeoStruct(variograms=v,nugget=0.5)``

        ``>>>pp_df = pyemu.pp_utils.pp_file_to_dataframe("hkpp.dat")``

        ``>>>ok = pyemu.utils.geostats.OrdinaryKrige(gs,pp_df)``

        ``>>>m = flopy.modflow.Modflow.load("mymodel.nam")``

        ``>>>df = ok.calc_factors_grid(m.sr,zone_array=m.bas6.ibound[0].array,``

        ``>>>                          var_filename="ok_var.dat")``

        ``>>>ok.to_grid_factor_file("factors.dat")``

def calc_factors_grid(self,spatial_reference,zone_array=None,minpts_interp=1,
                          maxpts_interp=20,search_radius=1.0e+10,verbose=False,
                          var_filename=None, forgive=False):
        """ calculate kriging factors (weights) for a structured grid.

        Parameters
        ----------
        spatial_reference : (flopy.utils.reference.SpatialReference)
            a spatial reference that describes the orientation and
            spatail projection of the the structured grid
        zone_array : (numpy.ndarray)
            an integer array of zones to use for kriging.  If not None,
            then point_data must also contain a "zone" column.  point_data
            entries with a zone value not found in zone_array will be skipped.
            If None, then all point_data will (potentially) be used for
            interpolating each grid node. Default is None
        minpts_interp : (int)
            minimum number of point_data entires to use for interpolation at
            a given grid node.  grid nodes with less than minpts_interp
            point_data found will be skipped (assigned np.NaN).  Defaut is 1
        maxpts_interp : (int)
            maximum number of point_data entries to use for interpolation at
            a given grid node.  A larger maxpts_interp will yield "smoother"
            interplation, but using a large maxpts_interp will slow the
            (already) slow kriging solution process and may lead to
            memory errors. Default is 20.
        search_radius : (float)
            the size of the region around a given grid node to search for
            point_data entries. Default is 1.0e+10
        verbose : (boolean)
            a flag to  echo process to stdout during the interpolatino process.
            Default is False
        var_filename : (str)
            a filename to save the kriging variance for each interpolated grid node.
            Default is None.
        forgive : (boolean)
            flag to continue if inversion of the kriging matrix failes at one or more
            grid nodes.  Inversion usually fails if the kriging matrix is singular,
            resulting from point_data entries closer than EPSILON distance.  If True,
            warnings are issued for each failed inversion.  If False, an exception
            is raised for failed matrix inversion.

        Returns
        -------
        df : pandas.DataFrame
            a dataframe with information summarizing the ordinary kriging
            process for each grid node

        Note
        ----
        this method calls OrdinaryKrige.calc_factors()


        Example
        -------
        ``>>>import flopy``

        ``>>>import pyemu``

        ``>>>v = pyemu.utils.geostats.ExpVario(a=1000,contribution=1.0)``

        ``>>>gs = pyemu.utils.geostats.GeoStruct(variograms=v,nugget=0.5)``

        ``>>>pp_df = pyemu.pp_utils.pp_file_to_dataframe("hkpp.dat")``

        ``>>>ok = pyemu.utils.geostats.OrdinaryKrige(gs,pp_df)``

        ``>>>m = flopy.modflow.Modflow.load("mymodel.nam")``

        ``>>>df = ok.calc_factors_grid(m.sr,zone_array=m.bas6.ibound[0].array,``

        ``>>>                          var_filename="ok_var.dat")``

        ``>>>ok.to_grid_factor_file("factors.dat")``

        """

        self.spatial_reference = spatial_reference
        self.interp_data = None
        #assert isinstance(spatial_reference,SpatialReference)
        try:
            x = self.spatial_reference.xcentergrid.copy()
            y = self.spatial_reference.ycentergrid.copy()
        except Exception as e:
            raise Exception("spatial_reference does not have proper attributes:{0}"\
                            .format(str(e)))

        if var_filename is not None:
                arr = np.zeros((self.spatial_reference.nrow,
                                self.spatial_reference.ncol)) - 1.0e+30

        # the simple case of no zone array: ignore point_data zones
        if zone_array is None:
            df = self.calc_factors(x.ravel(),y.ravel(),
                               minpts_interp=minpts_interp,
                               maxpts_interp=maxpts_interp,
                               search_radius=search_radius,
                               verbose=verbose, forgive=forgive)
            if var_filename is not None:
                arr = df.err_var.values.reshape(x.shape)
                np.savetxt(var_filename,arr,fmt="%15.6E")

        if zone_array is not None:
            assert zone_array.shape == x.shape
            if "zone" not in self.point_data.columns:
                warnings.warn("'zone' columns not in point_data, assigning generic zone",PyemuWarning)
                self.point_data.loc[:,"zone"] = 1
            pt_data_zones = self.point_data.zone.unique()
            dfs = []
            for pt_data_zone in pt_data_zones:
                if pt_data_zone not in zone_array:
                    warnings.warn("pt zone {0} not in zone array {1}, skipping".\
                                  format(pt_data_zone,np.unique(zone_array)),PyemuWarning)
                    continue
                xzone,yzone = x.copy(),y.copy()
                xzone[zone_array!=pt_data_zone] = np.NaN
                yzone[zone_array!=pt_data_zone] = np.NaN
                df = self.calc_factors(xzone.ravel(),yzone.ravel(),
                                       minpts_interp=minpts_interp,
                                       maxpts_interp=maxpts_interp,
                                       search_radius=search_radius,
                                       verbose=verbose,pt_zone=pt_data_zone,
                                       forgive=forgive)
                dfs.append(df)
                if var_filename is not None:
                    a = df.err_var.values.reshape(x.shape)
                    na_idx = np.isfinite(a)
                    arr[na_idx] = a[na_idx]
            if self.interp_data is None or self.interp_data.dropna().shape[0] == 0:
                raise Exception("no interpolation took place...something is wrong")
            df = pd.concat(dfs)
        if var_filename is not None:
            np.savetxt(var_filename,arr,fmt="%15.6E")
        return df