calculate ordinary kriging factors (weights) for the points
        represented by arguments x and y

        Parameters
        ----------
        x : (iterable of floats)
            x-coordinates to calculate kriging factors for
        y : (iterable of floats)
            y-coordinates to calculate kriging factors for
        minpts_interp : (int)
            minimum number of point_data entires to use for interpolation at
            a given x,y interplation point.  interpolation points 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 x,y interpolation point.  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 x,y interpolation point 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
        forgive : (boolean)
            flag to continue if inversion of the kriging matrix failes at one or more
            interpolation points.  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 interpolation points

def calc_factors(self,x,y,minpts_interp=1,maxpts_interp=20,
                     search_radius=1.0e+10,verbose=False,
                     pt_zone=None,forgive=False):
        """ calculate ordinary kriging factors (weights) for the points
        represented by arguments x and y

        Parameters
        ----------
        x : (iterable of floats)
            x-coordinates to calculate kriging factors for
        y : (iterable of floats)
            y-coordinates to calculate kriging factors for
        minpts_interp : (int)
            minimum number of point_data entires to use for interpolation at
            a given x,y interplation point.  interpolation points 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 x,y interpolation point.  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 x,y interpolation point 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
        forgive : (boolean)
            flag to continue if inversion of the kriging matrix failes at one or more
            interpolation points.  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 interpolation points


        """

        assert len(x) == len(y)

        # find the point data to use for each interp point
        sqradius = search_radius**2
        df = pd.DataFrame(data={'x':x,'y':y})
        inames,idist,ifacts,err_var = [],[],[],[]
        sill = self.geostruct.sill
        if pt_zone is None:
            ptx_array = self.point_data.x.values
            pty_array = self.point_data.y.values
            ptnames = self.point_data.name.values
        else:
            pt_data = self.point_data
            ptx_array = pt_data.loc[pt_data.zone==pt_zone,"x"].values
            pty_array = pt_data.loc[pt_data.zone==pt_zone,"y"].values
            ptnames = pt_data.loc[pt_data.zone==pt_zone,"name"].values
        #if verbose:
        print("starting interp point loop for {0} points".format(df.shape[0]))
        start_loop = datetime.now()
        for idx,(ix,iy) in enumerate(zip(df.x,df.y)):
            if np.isnan(ix) or np.isnan(iy): #if nans, skip
                inames.append([])
                idist.append([])
                ifacts.append([])
                err_var.append(np.NaN)
                continue
            if verbose:
                istart = datetime.now()
                print("processing interp point:{0} of {1}".format(idx,df.shape[0]))
            # if verbose == 2:
            #     start = datetime.now()
            #     print("calc ipoint dist...",end='')

            #  calc dist from this interp point to all point data...slow
            dist = pd.Series((ptx_array-ix)**2 + (pty_array-iy)**2,ptnames)
            dist.sort_values(inplace=True)
            dist = dist.loc[dist <= sqradius]

            # if too few points were found, skip
            if len(dist) < minpts_interp:
                inames.append([])
                idist.append([])
                ifacts.append([])
                err_var.append(sill)
                continue

            # only the maxpts_interp points
            dist = dist.iloc[:maxpts_interp].apply(np.sqrt)
            pt_names = dist.index.values
            # if one of the points is super close, just use it and skip
            if dist.min() <= EPSILON:
                ifacts.append([1.0])
                idist.append([EPSILON])
                inames.append([dist.idxmin()])
                err_var.append(self.geostruct.nugget)
                continue
            # if verbose == 2:
            #     td = (datetime.now()-start).total_seconds()
            #     print("...took {0}".format(td))
            #     start = datetime.now()
            #     print("extracting pt cov...",end='')

            #vextract the point-to-point covariance matrix
            point_cov = self.point_cov_df.loc[pt_names,pt_names]
            # if verbose == 2:
            #     td = (datetime.now()-start).total_seconds()
            #     print("...took {0}".format(td))
            #     print("forming ipt-to-point cov...",end='')

            # calc the interp point to points covariance
            interp_cov = self.geostruct.covariance_points(ix,iy,self.point_data.loc[pt_names,"x"],
                                                          self.point_data.loc[pt_names,"y"])

            if verbose == 2:
                td = (datetime.now()-start).total_seconds()
                print("...took {0} seconds".format(td))
                print("forming lin alg components...",end='')

            # form the linear algebra parts and solve
            d = len(pt_names) + 1 # +1 for lagrange mult
            A = np.ones((d,d))
            A[:-1,:-1] = point_cov.values
            A[-1,-1] = 0.0 #unbiaised constraint
            rhs = np.ones((d,1))
            rhs[:-1,0] = interp_cov
            # if verbose == 2:
            #     td = (datetime.now()-start).total_seconds()
            #     print("...took {0}".format(td))
            #     print("solving...",end='')
            # # solve
            try:
                facs = np.linalg.solve(A,rhs)
            except Exception as e:
                print("error solving for factors: {0}".format(str(e)))
                print("point:",ix,iy)
                print("dist:",dist)
                print("A:", A)
                print("rhs:", rhs)
                if forgive:
                    inames.append([])
                    idist.append([])
                    ifacts.append([])
                    err_var.append(np.NaN)
                    continue
                else:
                    raise Exception("error solving for factors:{0}".format(str(e)))
            assert len(facs) - 1 == len(dist)

            err_var.append(float(sill + facs[-1] - sum([f*c for f,c in zip(facs[:-1],interp_cov)])))
            inames.append(pt_names)

            idist.append(dist.values)
            ifacts.append(facs[:-1,0])
            # if verbose == 2:
            #     td = (datetime.now()-start).total_seconds()
            #     print("...took {0}".format(td))
            if verbose:
                td = (datetime.now()-istart).total_seconds()
                print("point took {0}".format(td))
        df["idist"] = idist
        df["inames"] = inames
        df["ifacts"] = ifacts
        df["err_var"] = err_var
        if pt_zone is None:
            self.interp_data = df
        else:
            if self.interp_data is None:
                self.interp_data = df
            else:
                self.interp_data = self.interp_data.append(df)
        td = (datetime.now() - start_loop).total_seconds()
        print("took {0} seconds".format(td))
        return df