def observe(self, obs_mesh, obs_V, output_type='r'):
        """
        Observes self on obs_mesh with observation variance obs_V.
        Output_type controls the information returned:

        'r' : returns information needed by Realization objects.
        'o' : returns information needed by function observe.
        's' : returns information needed by the Gaussian process
              submodel.
        """

        # Number of spatial dimensions.
        ndim = obs_mesh.shape[1]

        if self.ndim is not None:
            if not ndim == self.ndim:
                raise ValueError(
                    "Dimension of observation mesh is not equal to dimension of base mesh.")
        else:
            self.ndim = ndim

        # print ndim

        # =====================================
        # = If self hasn't been observed yet: =
        # =====================================
        if not self.observed:

            # If self has not been observed, get the Cholesky factor of self(obs_mesh, obs_mesh)
            # and the side information and store it.

            # Number of observation points so far is 0.
            N_old = 0
            N_new = obs_mesh.shape[0]

            if output_type == 's':
                U, C_eval = self.cholesky(
                    obs_mesh, nugget=obs_V, observed=False, return_eval_also=True)
                U_new = U
            else:
                U = self.cholesky(obs_mesh, nugget=obs_V, observed=False)

            # Upper-triangular Cholesky factor of self(obs_mesh, obs_mesh)
            self.full_Uo = U
            self.Uo = U

            # Pivots.
            piv_new = arange(N_new)
            self.full_piv = piv_new
            self.obs_piv = piv_new

            # Remember full observation mesh.
            self.full_obs_mesh = obs_mesh

            # relevant slice is the positive-definite indices, which get into
            # obs_mesh_*. See documentation.
            relevant_slice = self.obs_piv

            self.obs_mesh = obs_mesh
            self.obs_V = obs_V
            self.obs_len = N_new

        # =======================================
        # = If self has been observed already:  =
        # =======================================
        else:

            # If self has been observed, get the Cholesky factor of the _full_ observation mesh (new
            # and old observations) using continue_cholesky, along with side
            # information, and store it.

            # Number of observations so far.
            N_old = self.full_obs_mesh.shape[0]

            # Number of new observations.
            N_new = obs_mesh.shape[0]

            # Call to self.continue_cholesky.
            if output_type == 's':
                U, U_new, C_eval = self.continue_cholesky(x=obs_mesh,
                                                          x_old=self.full_obs_mesh,
                                                          U_old=self.full_Uo,
                                                          observed=False,
                                                          nugget=obs_V,
                                                          return_eval_also=True)

            else:
                U = self.continue_cholesky(x=obs_mesh,
                                           x_old=self.full_obs_mesh,
                                           U_old=self.full_Uo,
                                           observed=False,
                                           nugget=obs_V)

            # Full Cholesky factor of self(obs_mesh, obs_mesh), where obs_mesh
            # is the combined observation mesh.
            self.full_Uo = U

            # Square upper-triangular Cholesky factor of self(obs_mesh_*,
            # obs_mesh_*). See documentation.
            self.Uo = self.full_Uo

            # Pivots.
            piv_new = arange(N_old + N_new)
            self.obs_piv = piv_new
            self.full_piv = piv_new

            # Concatenate old and new observation meshes.
            self.full_obs_mesh = vstack((self.full_obs_mesh, obs_mesh))
            self.obs_mesh = self.full_obs_mesh
            self.obs_V = hstack((self.obs_V, obs_V))

            # Length of obs_mesh_*.
            self.obs_len = N_old + N_new

        self.observed = True
        # Output expected by Realization
        if output_type == 'r':
            return slice(None, None, None), obs_mesh, self.full_Uo[
                N_old:N_new + N_old, N_old:N_new + N_old], self.full_Uo[:N_old, N_old:N_new + N_old]

        # Ouptut expected by observe
        if output_type == 'o':
            return slice(None, None, None), obs_mesh

        # Output expected by the GP submodel
        if output_type == 's':
            return U_new, C_eval, self.full_Uo[:N_old, N_old:N_new + N_old]
Observes self on obs_mesh with observation variance obs_V.
        Output_type controls the information returned:

        'r' : returns information needed by Realization objects.
        'o' : returns information needed by function observe.
        's' : returns information needed by the Gaussian process
              submodel.