diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py
index 2eb980c4a..6bb0c78ce 100644
--- a/gstools/covmodel/base.py
+++ b/gstools/covmodel/base.py
@@ -39,7 +39,7 @@
 
 
 class CovModel(six.with_metaclass(InitSubclassMeta)):
-    """Base class for the GSTools covariance models
+    """Base class for the GSTools covariance models.
 
     Parameters
     ----------
@@ -184,7 +184,7 @@ def __init__(
     ###########################################################################
 
     def __init_subclass__(cls):
-        r"""Initialize gstools covariance model
+        r"""Initialize gstools covariance model.
 
         Warnings
         --------
@@ -206,8 +206,9 @@ def __init_subclass__(cls):
         Best practice is to use the ``correlation`` function!
         """
         # overrid one of these ################################################
+
         def variogram(self, r):
-            r"""Isotropic variogram of the model
+            r"""Isotropic variogram of the model.
 
             Given by: :math:`\gamma\left(r\right)=
             \sigma^2\cdot\left(1-\mathrm{cor}\left(r\right)\right)+n`
@@ -217,7 +218,7 @@ def variogram(self, r):
             return self.var - self.covariance(r) + self.nugget
 
         def covariance(self, r):
-            r"""Covariance of the model
+            r"""Covariance of the model.
 
             Given by: :math:`C\left(r\right)=
             \sigma^2\cdot\mathrm{cor}\left(r\right)`
@@ -227,7 +228,7 @@ def covariance(self, r):
             return self.var * self.correlation(r)
 
         def correlation(self, r):
-            r"""correlation function (or normalized covariance) of the model
+            r"""Correlation function (or normalized covariance) of the model.
 
             Given by: :math:`\mathrm{cor}\left(r\right)`
 
@@ -237,7 +238,7 @@ def correlation(self, r):
             return 1.0 - self.variogram_normed(r)
 
         def variogram_normed(self, r):
-            r"""Normalized variogram of the model
+            r"""Normalized-variogram of the model.
 
             Given by: :math:`\tilde{\gamma}\left(r\right)=
             1-\mathrm{cor}\left(r\right)`
@@ -300,31 +301,31 @@ def variogram_normed(self, r):
     ### special variogram functions ###########################################
 
     def cov_nugget(self, r):
-        r"""Covariance of the model respecting the nugget at r=0
+        r"""Covariance of the model respecting the nugget at r=0.
 
         Given by: :math:`C\left(r\right)=
         \sigma^2\cdot\mathrm{cor}\left(r\right)`
 
         Where :math:`\mathrm{cor}(r)` is the correlation function.
         """
-        r = np.array(np.abs(r), dtype=float)
+        r = np.array(np.abs(r), dtype=np.double)
         r_gz = np.logical_not(np.isclose(r, 0))
-        res = np.empty_like(r, dtype=float)
+        res = np.empty_like(r, dtype=np.double)
         res[r_gz] = self.covariance(r[r_gz])
         res[np.logical_not(r_gz)] = self.sill
         return res
 
     def vario_nugget(self, r):
-        r"""Isotropic variogram of the model respecting the nugget at r=0
+        r"""Isotropic variogram of the model respecting the nugget at r=0.
 
         Given by: :math:`\gamma\left(r\right)=
         \sigma^2\cdot\left(1-\mathrm{cor}\left(r\right)\right)+n`
 
         Where :math:`\mathrm{cor}(r)` is the correlation function.
         """
-        r = np.array(np.abs(r), dtype=float)
+        r = np.array(np.abs(r), dtype=np.double)
         r_gz = np.logical_not(np.isclose(r, 0))
-        res = np.empty_like(r, dtype=float)
+        res = np.empty_like(r, dtype=np.double)
         res[r_gz] = self.variogram(r[r_gz])
         res[np.logical_not(r_gz)] = 0.0
         return res
@@ -334,7 +335,7 @@ def vario_nugget(self, r):
     ###########################################################################
 
     def pykrige_vario(self, args=None, r=0):
-        r"""Isotropic variogram of the model for pykrige
+        r"""Isotropic variogram of the model for pykrige.
 
         Given by: :math:`\gamma\left(r\right)=
         \sigma^2\cdot\left(1-\mathrm{cor}\left(r\right)\right)+n`
@@ -345,56 +346,56 @@ def pykrige_vario(self, args=None, r=0):
 
     @property
     def pykrige_anis(self):
-        """2D anisotropy ratio for pykrige"""
+        """2D anisotropy ratio for pykrige."""
         if self.dim == 2:
             return 1 / self.anis[0]
         return 1.0
 
     @property
     def pykrige_anis_y(self):
-        """3D anisotropy ratio in y direction for pykrige"""
+        """3D anisotropy ratio in y direction for pykrige."""
         if self.dim >= 2:
             return 1 / self.anis[0]
         return 1.0
 
     @property
     def pykrige_anis_z(self):
-        """3D anisotropy ratio in z direction for pykrige"""
+        """3D anisotropy ratio in z direction for pykrige."""
         if self.dim == 3:
             return 1 / self.anis[1]
         return 1.0
 
     @property
     def pykrige_angle(self):
-        """2D rotation angle for pykrige"""
+        """2D rotation angle for pykrige."""
         if self.dim == 2:
             return -self.angles[0] / np.pi * 180
         return 0.0
 
     @property
     def pykrige_angle_z(self):
-        """3D rotation angle around z for pykrige"""
+        """3D rotation angle around z for pykrige."""
         if self.dim >= 2:
             return -self.angles[0] / np.pi * 180
         return 0.0
 
     @property
     def pykrige_angle_y(self):
-        """3D rotation angle around y for pykrige"""
+        """3D rotation angle around y for pykrige."""
         if self.dim == 3:
             return -self.angles[1] / np.pi * 180
         return 0.0
 
     @property
     def pykrige_angle_x(self):
-        """3D rotation angle around x for pykrige"""
+        """3D rotation angle around x for pykrige."""
         if self.dim == 3:
             return -self.angles[2] / np.pi * 180
         return 0.0
 
     @property
     def pykrige_kwargs(self):
-        """keyword arguments for pykrige routines"""
+        """Keyword arguments for pykrige routines."""
         kwargs = {
             "variogram_model": "custom",
             "variogram_parameters": [],
@@ -423,20 +424,21 @@ def pykrige_kwargs(self):
     ###########################################################################
 
     def default_opt_arg(self):
-        """Here you can provide a dictionary with default values for
-        the optional arguments."""
+        """Provide default optional arguments by the user.
+
+        Should be given as a dictionary.
+        """
         return {}
 
     def default_opt_arg_bounds(self):
-        """Here you can provide a dictionary with default boundaries for
-        the optional arguments."""
+        """Provide default boundaries for optional arguments."""
         res = {}
         for opt in self.opt_arg:
             res[opt] = [0.0, 1000.0]
         return res
 
     def check_opt_arg(self):
-        """Here you can run checks for the optional arguments
+        """Run checks for the optional arguments.
 
         This is in addition to the bound-checks
 
@@ -449,23 +451,26 @@ def check_opt_arg(self):
         pass
 
     def fix_dim(self):
-        """Set a fix dimension for the model"""
+        """Set a fix dimension for the model."""
         return None
 
     def var_factor(self):
-        """Optional factor for the variance"""
+        """Factor for the variance."""
         return 1.0
 
     ### calculation of different scales #######################################
 
     def _calc_integral_scale(self):
-        """calculate the integral scale of the isotrope model"""
+        """Calculate the integral scale of the isotrope model."""
         self._integral_scale = integral(self.correlation, 0, np.inf)[0]
         return self._integral_scale
 
     def percentile_scale(self, per=0.9):
-        """calculate the distance, where the given percentile of the variance
-        is reached by the variogram"""
+        """Calculate the percentile scale of the isotrope model.
+
+        This is the distance, where the given percentile of the variance
+        is reached by the variogram
+        """
         # check the given percentile
         if not 0.0 < per < 1.0:
             raise ValueError(
@@ -485,7 +490,7 @@ def curve(x):
 
     def spectrum(self, k):
         r"""
-        The spectrum of the covariance model.
+        Spectrum of the covariance model.
 
         This is given by:
 
@@ -505,12 +510,12 @@ def spectrum(self, k):
         k : :class:`float`
             Radius of the phase: :math:`k=\left\Vert\mathbf{k}\right\Vert`
         """
-        k = np.array(np.abs(k), dtype=float)
+        k = np.array(np.abs(k), dtype=np.double)
         if self.dim > 1:
             res = self._sft.transform(self.covariance, k, ret_err=False)
         else:
             k_gz = np.logical_not(np.isclose(k, 0))
-            res = np.empty_like(k, dtype=float)
+            res = np.empty_like(k, dtype=np.double)
             res[k_gz] = self._sft.transform(
                 self.covariance, k[k_gz], ret_err=False
             )
@@ -527,7 +532,7 @@ def spectrum(self, k):
 
     def spectral_density(self, k):
         r"""
-        The spectral density of the covariance model.
+        Spectral density of the covariance model.
 
         This is given by:
 
@@ -543,14 +548,12 @@ def spectral_density(self, k):
         return self.spectrum(k) / self.var
 
     def spectral_rad_pdf(self, r):
-        """
-        The radial spectral density of the model depending on the dimension
-        """
-        r = np.array(np.abs(r), dtype=float)
+        """Radial spectral density of the model."""
+        r = np.array(np.abs(r), dtype=np.double)
         if self.dim > 1:
             r_gz = np.logical_not(np.isclose(r, 0))
             # to prevent numerical errors, we just calculate where r>0
-            res = np.zeros_like(r, dtype=float)
+            res = np.zeros_like(r, dtype=np.double)
             res[r_gz] = rad_fac(self.dim, r[r_gz]) * self.spectral_density(
                 r[r_gz]
             )
@@ -563,28 +566,26 @@ def spectral_rad_pdf(self, r):
         return res
 
     def ln_spectral_rad_pdf(self, r):
-        """
-        The log radial spectral density of the model depending on the dimension
-        """
+        """Log radial spectral density of the model."""
         spec = np.array(self.spectral_rad_pdf(r))
         spec_gz = np.logical_not(np.isclose(spec, 0))
-        res = np.full_like(spec, -np.inf, dtype=float)
+        res = np.full_like(spec, -np.inf, dtype=np.double)
         res[spec_gz] = np.log(spec[spec_gz])
         return res
 
     def _has_cdf(self):
-        """State if a cdf is defined with 'spectral_rad_cdf'"""
+        """State if a cdf is defined with 'spectral_rad_cdf'."""
         return hasattr(self, "spectral_rad_cdf")
 
     def _has_ppf(self):
-        """State if a ppf is defined with 'spectral_rad_ppf'"""
+        """State if a ppf is defined with 'spectral_rad_ppf'."""
         return hasattr(self, "spectral_rad_ppf")
 
     ### fitting routine #######################################################
 
     def fit_variogram(self, x_data, y_data, **para_deselect):
         """
-        fit the isotropic variogram-model to given data
+        Fiting the isotropic variogram-model to given data.
 
         Parameters
         ----------
@@ -621,7 +622,7 @@ def fit_variogram(self, x_data, y_data, **para_deselect):
 
         # we need arg1, otherwise curve_fit throws an error (bug?!)
         def curve(x, arg1, *args):
-            """dummy function for the variogram"""
+            """Adapted Variogram function."""
             args = (arg1,) + args
             para_skip = 0
             opt_skip = 0
@@ -699,8 +700,10 @@ def curve(x, arg1, *args):
     ### bounds setting and checks #############################################
 
     def default_arg_bounds(self):
-        """Here you can provide a dictionary with default boundaries for
-        the standard arguments."""
+        """Provide default boundaries for arguments.
+
+        Given as a dictionary.
+        """
         res = {
             "var": (0.0, 100.0, "oc"),
             "len_scale": (0.0, 1000.0, "oo"),
@@ -709,7 +712,7 @@ def default_arg_bounds(self):
         return res
 
     def set_arg_bounds(self, **kwargs):
-        r"""Set bounds for the parameters of the model
+        r"""Set bounds for the parameters of the model.
 
         Parameters
         ----------
@@ -741,7 +744,7 @@ def set_arg_bounds(self, **kwargs):
                 self.nugget_bounds = kwargs[opt]
 
     def check_arg_bounds(self):
-        """Here the arguments are checked to be within the given bounds"""
+        """Check arguments to be within the given bounds."""
         # check var, len_scale, nugget and optional-arguments
         for arg in self.arg_bounds:
             bnd = list(self.arg_bounds[arg])
@@ -789,7 +792,7 @@ def check_arg_bounds(self):
 
     @property
     def var_bounds(self):
-        """:class:`list`: Bounds for the variance
+        """:class:`list`: Bounds for the variance.
 
         Notes
         -----
@@ -815,7 +818,7 @@ def var_bounds(self, bounds):
 
     @property
     def len_scale_bounds(self):
-        """:class:`list`: Bounds for the lenght scale
+        """:class:`list`: Bounds for the lenght scale.
 
         Notes
         -----
@@ -841,7 +844,7 @@ def len_scale_bounds(self, bounds):
 
     @property
     def nugget_bounds(self):
-        """:class:`list`: Bounds for the nugget
+        """:class:`list`: Bounds for the nugget.
 
         Notes
         -----
@@ -867,7 +870,7 @@ def nugget_bounds(self, bounds):
 
     @property
     def opt_arg_bounds(self):
-        """:class:`dict`: Bounds for the optional arguments
+        """:class:`dict`: Bounds for the optional arguments.
 
         Notes
         -----
@@ -883,7 +886,7 @@ def opt_arg_bounds(self):
 
     @property
     def arg_bounds(self):
-        """:class:`dict`: Bounds for all parameters
+        """:class:`dict`: Bounds for all parameters.
 
         Notes
         -----
@@ -942,9 +945,10 @@ def var(self, var):
 
     @property
     def var_raw(self):
-        """:class:`float`: The raw variance of the model without factor
+        """:class:`float`: The raw variance of the model without factor.
 
-        (See. CovModel.var_factor)"""
+        (See. CovModel.var_factor)
+        """
         return self._var
 
     @var_raw.setter
@@ -988,7 +992,7 @@ def anis(self, anis):
 
     @property
     def angles(self):
-        """:class:`numpy.ndarray`: The rotation angles (in rad) of the model."""
+        """:class:`numpy.ndarray`: Rotation angles (in rad) of the model."""
         return self._angles
 
     @angles.setter
@@ -1027,8 +1031,7 @@ def integral_scale(self, integral_scale):
 
     @property
     def hankel_kw(self):
-        """:class:`dict`: Keywords for :any:`hankel.SymmetricFourierTransform`
-        """
+        """:class:`dict`: :any:`hankel.SymmetricFourierTransform` kwargs."""
         return self._hankel_kw
 
     @hankel_kw.setter
@@ -1041,8 +1044,10 @@ def hankel_kw(self, hankel_kw):
 
     @property
     def dist_func(self):
-        """:class:`tuple` of :any:`callable`:
-        pdf, cdf and ppf from the radial spectral density"""
+        """:class:`tuple` of :any:`callable`: pdf, cdf and ppf.
+
+        Spectral distribution info from the model.
+        """
         pdf = self.spectral_rad_pdf
         cdf = None
         ppf = None
@@ -1054,12 +1059,12 @@ def dist_func(self):
 
     @property
     def has_cdf(self):
-        """:class:`bool`: State if a cdf is defined by the user"""
+        """:class:`bool`: State if a cdf is defined by the user."""
         return self._has_cdf()
 
     @property
     def has_ppf(self):
-        """:class:`bool`: State if a ppf is defined by the user"""
+        """:class:`bool`: State if a ppf is defined by the user."""
         return self._has_ppf()
 
     @property
@@ -1075,12 +1080,12 @@ def sill(self):
 
     @property
     def arg(self):
-        """:class:`list` of :class:`str`: Names of all arguments"""
+        """:class:`list` of :class:`str`: Names of all arguments."""
         return ["var", "len_scale", "nugget", "anis", "angles"] + self._opt_arg
 
     @property
     def opt_arg(self):
-        """:class:`list` of :class:`str`: Names of the optional arguments"""
+        """:class:`list` of :class:`str`: Names of the optional arguments."""
         return self._opt_arg
 
     @property
@@ -1094,7 +1099,7 @@ def len_scale_vec(self):
             * ``len_scale_vec[1] = len_scale*anis[0]``
             * ``len_scale_vec[2] = len_scale*anis[1]``
         """
-        res = np.zeros(self.dim, dtype=float)
+        res = np.zeros(self.dim, dtype=np.double)
         res[0] = self.len_scale
         for i in range(1, self._dim):
             res[i] = self.len_scale * self.anis[i - 1]
@@ -1111,7 +1116,7 @@ def integral_scale_vec(self):
             * ``integral_scale_vec[1] = integral_scale*anis[0]``
             * ``integral_scale_vec[2] = integral_scale*anis[1]``
         """
-        res = np.zeros(self.dim, dtype=float)
+        res = np.zeros(self.dim, dtype=np.double)
         res[0] = self.integral_scale
         for i in range(1, self.dim):
             res[i] = self.integral_scale * self.anis[i - 1]
@@ -1119,12 +1124,13 @@ def integral_scale_vec(self):
 
     @property
     def name(self):
-        """:class:`str`: The name of the CovModel class"""
+        """:class:`str`: The name of the CovModel class."""
         return self.__class__.__name__
 
     ### magic methods #########################################################
 
     def __eq__(self, other):
+        """Compare CovModels."""
         if not isinstance(other, CovModel):
             return False
         # prevent attribute error in opt_arg if the are not equal
@@ -1146,12 +1152,15 @@ def __eq__(self, other):
         return equal
 
     def __ne__(self, other):
+        """Compare CovModels."""
         return not self.__eq__(other)
 
     def __str__(self):
+        """Return String representation."""
         return self.__repr__()
 
     def __repr__(self):
+        """Return String representation."""
         opt_str = ""
         for opt in self.opt_arg:
             opt_str += ", " + opt + "={}".format(getattr(self, opt))
diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py
index cd247db73..3bdc9bfc1 100644
--- a/gstools/covmodel/models.py
+++ b/gstools/covmodel/models.py
@@ -42,7 +42,7 @@
 
 
 class Gaussian(CovModel):
-    r"""The Gaussian covariance model
+    r"""The Gaussian covariance model.
 
     Notes
     -----
@@ -54,13 +54,13 @@ class Gaussian(CovModel):
     """
 
     def correlation(self, r):
-        r"""Gaussian correlation function
+        r"""Gaussian correlation function.
 
         .. math::
            \mathrm{cor}(r) =
-           \exp\left(- \frac{\pi}{4} \cdot \left(\frac{r}{\ell}\right)^2\right)
+           \exp\left(- \frac{\pi}{4}\cdot \left(\frac{r}{\ell}\right)^2\right)
        """
-        r = np.array(np.abs(r), dtype=float)
+        r = np.array(np.abs(r), dtype=np.double)
         return np.exp(-np.pi / 4 * (r / self.len_scale) ** 2)
 
     def spectrum(self, k):
@@ -71,7 +71,7 @@ def spectrum(self, k):
         )
 
     def spectral_rad_cdf(self, r):
-        """The cdf of the radial spectral density"""
+        """Radial spectral cdf."""
         if self.dim == 1:
             return sps.erf(self.len_scale * r / np.sqrt(np.pi))
         if self.dim == 2:
@@ -85,7 +85,7 @@ def spectral_rad_cdf(self, r):
         return None
 
     def spectral_rad_ppf(self, u):
-        """The ppf of the radial spectral density"""
+        """Radial spectral ppf."""
         if self.dim == 1:
             return sps.erfinv(u) * np.sqrt(np.pi) / self.len_scale
         if self.dim == 2:
@@ -93,14 +93,17 @@ def spectral_rad_ppf(self, u):
         return None
 
     def _has_ppf(self):
-        """ppf for 3 dimensions is not analytical given"""
+        """State if ppf is given.
+
+        ppf for 3 dimensions is not analytical
+        """
         # since the ppf is not analytical for dim=3, we have to state that
         if self.dim == 3:
             return False
         return True
 
     def calc_integral_scale(self):
-        """The integral scale of the gaussian model is the length scale"""
+        """Integral scale of the gaussian model is the length scale."""
         return self.len_scale
 
 
@@ -108,7 +111,7 @@ def calc_integral_scale(self):
 
 
 class Exponential(CovModel):
-    r"""The Exponential covariance model
+    r"""The Exponential covariance model.
 
     Notes
     -----
@@ -120,13 +123,13 @@ class Exponential(CovModel):
     """
 
     def correlation(self, r):
-        r"""Exponential correlation function
+        r"""Exponential correlation function.
 
         .. math::
            \mathrm{cor}(r) =
            \exp\left(- \frac{r}{\ell} \right)
        """
-        r = np.array(np.abs(r), dtype=float)
+        r = np.array(np.abs(r), dtype=np.double)
         return np.exp(-1 * r / self.len_scale)
 
     def spectrum(self, k):
@@ -139,7 +142,7 @@ def spectrum(self, k):
         )
 
     def spectral_rad_cdf(self, r):
-        """The cdf of the radial spectral density"""
+        """Radial spectral cdf."""
         if self.dim == 1:
             return np.arctan(r * self.len_scale) * 2 / np.pi
         if self.dim == 2:
@@ -156,7 +159,7 @@ def spectral_rad_cdf(self, r):
         return None
 
     def spectral_rad_ppf(self, u):
-        """The ppf of the radial spectral density"""
+        """Radial spectral ppf."""
         if self.dim == 1:
             return np.tan(np.pi / 2 * u) / self.len_scale
         if self.dim == 2:
@@ -164,14 +167,17 @@ def spectral_rad_ppf(self, u):
         return None
 
     def _has_ppf(self):
-        """ppf for 3 dimensions is not analytical"""
+        """State if ppf is given.
+
+        ppf for 3 dimensions is not analytical
+        """
         # since the ppf is not analytical for dim=3, we have to state that
         if self.dim == 3:
             return False
         return True
 
     def calc_integral_scale(self):
-        """The integral scale of the exponential model is the length scale"""
+        """Integral scale of the exponential model is the length scale."""
         return self.len_scale
 
 
@@ -179,7 +185,7 @@ def calc_integral_scale(self):
 
 
 class Spherical(CovModel):
-    r"""The Spherical covariance model
+    r"""The Spherical covariance model.
 
     Notes
     -----
@@ -196,7 +202,7 @@ class Spherical(CovModel):
     """
 
     def correlation(self, r):
-        r"""Spherical correlation function
+        r"""Spherical correlation function.
 
         .. math::
            \mathrm{cor}(r) =
@@ -207,7 +213,7 @@ def correlation(self, r):
            0 & r\geq\ell
            \end{cases}
         """
-        r = np.array(np.abs(r), dtype=float)
+        r = np.array(np.abs(r), dtype=np.double)
         res = np.zeros_like(r)
         res[r < self.len_scale] = (
             1.0
@@ -218,7 +224,7 @@ def correlation(self, r):
 
 
 class SphericalRescal(CovModel):
-    r"""The rescaled Spherical covariance model
+    r"""The rescaled Spherical covariance model.
 
     Notes
     -----
@@ -235,7 +241,7 @@ class SphericalRescal(CovModel):
     """
 
     def correlation(self, r):
-        r"""Rescaled Spherical correlation function
+        r"""Rescaled Spherical correlation function.
 
         .. math::
            \mathrm{cor}(r) =
@@ -246,7 +252,7 @@ def correlation(self, r):
            0 & r\geq\frac{8}{3}\ell
            \end{cases}
         """
-        r = np.array(np.abs(r), dtype=float)
+        r = np.array(np.abs(r), dtype=np.double)
         res = np.zeros_like(r)
         res[r < 8 / 3 * self.len_scale] = (
             1.0
@@ -258,7 +264,7 @@ def correlation(self, r):
         return res
 
     def calc_integral_scale(self):
-        """The integral scale of this spherical model is the length scale"""
+        """Integral scale of this spherical model is the length scale."""
         return self.len_scale
 
 
@@ -266,7 +272,7 @@ def calc_integral_scale(self):
 
 
 class Rational(CovModel):
-    r"""The rational quadratic covariance model
+    r"""The rational quadratic covariance model.
 
     Notes
     -----
@@ -289,7 +295,7 @@ class Rational(CovModel):
     """
 
     def default_opt_arg(self):
-        """The defaults for the optional arguments:
+        """Defaults for the optional arguments.
 
             * ``{"alpha": 1.0}``
 
@@ -301,7 +307,7 @@ def default_opt_arg(self):
         return {"alpha": 1.0}
 
     def default_opt_arg_bounds(self):
-        """The defaults boundaries for the optional arguments:
+        """Defaults for boundaries of the optional arguments.
 
             * ``{"alpha": [0.5, inf]}``
 
@@ -313,14 +319,14 @@ def default_opt_arg_bounds(self):
         return {"alpha": [0.5, np.inf]}
 
     def correlation(self, r):
-        r"""Rational correlation function
+        r"""Rational correlation function.
 
         .. math::
            \mathrm{cor}(r) =
            \left(1 + \frac{1}{2\alpha} \cdot
            \left(\frac{r}{\ell}\right)^2\right)^{-\alpha}
         """
-        r = np.array(np.abs(r), dtype=float)
+        r = np.array(np.abs(r), dtype=np.double)
         return np.power(
             1 + 0.5 / self.alpha * (r / self.len_scale) ** 2, -self.alpha
         )
@@ -330,7 +336,7 @@ def correlation(self, r):
 
 
 class Stable(CovModel):
-    r"""The stable covariance model
+    r"""The stable covariance model.
 
     Notes
     -----
@@ -352,7 +358,7 @@ class Stable(CovModel):
     """
 
     def default_opt_arg(self):
-        """The defaults for the optional arguments:
+        """Defaults for the optional arguments.
 
             * ``{"alpha": 1.5}``
 
@@ -364,7 +370,7 @@ def default_opt_arg(self):
         return {"alpha": 1.5}
 
     def default_opt_arg_bounds(self):
-        """The defaults boundaries for the optional arguments:
+        """Defaults for boundaries of the optional arguments.
 
             * ``{"alpha": [0, 2, "oc"]}``
 
@@ -376,7 +382,7 @@ def default_opt_arg_bounds(self):
         return {"alpha": [0, 2, "oc"]}
 
     def check_opt_arg(self):
-        """Checks for the optional arguments
+        """Check the optional arguments.
 
         Warns
         -----
@@ -386,18 +392,18 @@ def check_opt_arg(self):
         """
         if self.alpha < 0.3:
             warnings.warn(
-                "TPLStable: parameter 'alpha' is < 0.3, "
+                "Stable: parameter 'alpha' is < 0.3, "
                 + "count with unstable results"
             )
 
     def correlation(self, r):
-        r"""Stable correlation function
+        r"""Stable correlation function.
 
         .. math::
            \mathrm{cor}(r) =
            \exp\left(- \left(\frac{r}{\ell}\right)^{\alpha}\right)
         """
-        r = np.array(np.abs(r), dtype=float)
+        r = np.array(np.abs(r), dtype=np.double)
         return np.exp(-np.power(r / self.len_scale, self.alpha))
 
 
@@ -405,7 +411,7 @@ def correlation(self, r):
 
 
 class Matern(CovModel):
-    r"""The Matérn covariance model
+    r"""The Matérn covariance model.
 
     Notes
     -----
@@ -432,7 +438,7 @@ class Matern(CovModel):
     """
 
     def default_opt_arg(self):
-        """The defaults for the optional arguments:
+        """Defaults for the optional arguments.
 
             * ``{"nu": 1.0}``
 
@@ -444,7 +450,7 @@ def default_opt_arg(self):
         return {"nu": 1.0}
 
     def default_opt_arg_bounds(self):
-        """The defaults boundaries for the optional arguments:
+        """Defaults for boundaries of the optional arguments.
 
             * ``{"nu": [0.5, 60.0, "cc"]}``
 
@@ -456,7 +462,7 @@ def default_opt_arg_bounds(self):
         return {"nu": [0.5, 60.0, "cc"]}
 
     def check_opt_arg(self):
-        """Checks for the optional arguments
+        """Check for the optional arguments.
 
         Warns
         -----
@@ -465,12 +471,12 @@ def check_opt_arg(self):
         """
         if self.nu > 50.0:
             warnings.warn(
-                "Mat: parameter 'nu' is > 50, "
+                "Matern: parameter 'nu' is > 50, "
                 + "calculations most likely get unstable here"
             )
 
     def correlation(self, r):
-        r"""Matérn correlation function
+        r"""Matérn correlation function.
 
         .. math::
            \mathrm{cor}(r) =
@@ -478,7 +484,7 @@ def correlation(self, r):
            \left(\sqrt{2\nu}\cdot\frac{r}{\ell}\right)^{\nu} \cdot
            \mathrm{K}_{\nu}\left(\sqrt{2\nu}\cdot\frac{r}{\ell}\right)
         """
-        r = np.array(np.abs(r), dtype=float)
+        r = np.array(np.abs(r), dtype=np.double)
         r_gz = r[r > 0.0]
         res = np.ones_like(r)
         with np.errstate(over="ignore", invalid="ignore"):
@@ -500,7 +506,7 @@ def correlation(self, r):
 
 
 class MaternRescal(CovModel):
-    r"""The rescaled Matérn covariance model
+    r"""The rescaled Matérn covariance model.
 
     Notes
     -----
@@ -530,7 +536,7 @@ class MaternRescal(CovModel):
     """
 
     def default_opt_arg(self):
-        """The defaults for the optional arguments:
+        """Defaults for the optional arguments.
 
             * ``{"nu": 1.0}``
 
@@ -542,7 +548,7 @@ def default_opt_arg(self):
         return {"nu": 1.0}
 
     def default_opt_arg_bounds(self):
-        """The defaults boundaries for the optional arguments:
+        """Defaults for boundaries of the optional arguments.
 
             * ``{"nu": [0.5, 60.0, "cc"]}``
 
@@ -554,7 +560,7 @@ def default_opt_arg_bounds(self):
         return {"nu": [0.5, 60.0, "cc"]}
 
     def check_opt_arg(self):
-        """Checks for the optional arguments
+        """Check for the optional arguments.
 
         Warns
         -----
@@ -563,12 +569,12 @@ def check_opt_arg(self):
         """
         if self.nu > 50.0:
             warnings.warn(
-                "Mat: parameter 'nu' is > 50, "
+                "MaternRescale: parameter 'nu' is > 50, "
                 + "calculations most likely get unstable here"
             )
 
     def correlation(self, r):
-        r"""Rescaled Matérn correlation function
+        r"""Rescaled Matérn correlation function.
 
         .. math::
            \mathrm{cor}(r) =
@@ -578,7 +584,7 @@ def correlation(self, r):
            \mathrm{K}_{\nu}\left(\frac{\pi}{B\left(\nu,\frac{1}{2}\right)}
            \cdot\frac{r}{\ell}\right)
         """
-        r = np.array(np.abs(r), dtype=float)
+        r = np.array(np.abs(r), dtype=np.double)
         r_gz = r[r > 0.0]
         res = np.ones_like(r)
         with np.errstate(over="ignore", invalid="ignore"):
@@ -605,7 +611,7 @@ def correlation(self, r):
 
 
 class Linear(CovModel):
-    r"""The bounded linear covariance model
+    r"""The bounded linear covariance model.
 
     Notes
     -----
@@ -621,7 +627,7 @@ class Linear(CovModel):
     """
 
     def correlation(self, r):
-        r"""Linear correlation function
+        r"""Linear correlation function.
 
         .. math::
            \mathrm{cor}(r) =
@@ -631,7 +637,7 @@ def correlation(self, r):
            0 & r\geq\ell
            \end{cases}
         """
-        r = np.array(np.abs(r), dtype=float)
+        r = np.array(np.abs(r), dtype=np.double)
         res = np.zeros_like(r)
         res[r < self.len_scale] = 1.0 - r[r < self.len_scale] / self.len_scale
         return res
diff --git a/gstools/covmodel/plot.py b/gstools/covmodel/plot.py
index aa985819e..1d5cc7329 100644
--- a/gstools/covmodel/plot.py
+++ b/gstools/covmodel/plot.py
@@ -34,7 +34,7 @@
 
 
 def plot_variogram(model, x_min=0.0, x_max=None):
-    """plot variogram of a given CovModel"""
+    """Plot variogram of a given CovModel."""
     if x_max is None:
         x_max = 3 * model.integral_scale
     x_s = np.linspace(x_min, x_max)
@@ -44,7 +44,7 @@ def plot_variogram(model, x_min=0.0, x_max=None):
 
 
 def plot_covariance(model, x_min=0.0, x_max=None):
-    """plot covariance of a given CovModel"""
+    """Plot covariance of a given CovModel."""
     if x_max is None:
         x_max = 3 * model.integral_scale
     x_s = np.linspace(x_min, x_max)
@@ -54,7 +54,7 @@ def plot_covariance(model, x_min=0.0, x_max=None):
 
 
 def plot_correlation(model, x_min=0.0, x_max=None):
-    """plot correlation function of a given CovModel"""
+    """Plot correlation function of a given CovModel."""
     if x_max is None:
         x_max = 3 * model.integral_scale
     x_s = np.linspace(x_min, x_max)
@@ -64,7 +64,7 @@ def plot_correlation(model, x_min=0.0, x_max=None):
 
 
 def plot_variogram_normed(model, x_min=0.0, x_max=None):
-    """plot normalized variogram of a given CovModel"""
+    """Plot normalized variogram of a given CovModel."""
     if x_max is None:
         x_max = 3 * model.integral_scale
     x_s = np.linspace(x_min, x_max)
@@ -76,7 +76,7 @@ def plot_variogram_normed(model, x_min=0.0, x_max=None):
 
 
 def plot_spectrum(model, x_min=0.0, x_max=None):
-    """plot specturm of a given CovModel"""
+    """Plot specturm of a given CovModel."""
     if x_max is None:
         x_max = 3 / model.integral_scale
     x_s = np.linspace(x_min, x_max)
@@ -90,7 +90,7 @@ def plot_spectrum(model, x_min=0.0, x_max=None):
 
 
 def plot_spectral_density(model, x_min=0.0, x_max=None):
-    """plot spectral density of a given CovModel"""
+    """Plot spectral density of a given CovModel."""
     if x_max is None:
         x_max = 3 / model.integral_scale
     x_s = np.linspace(x_min, x_max)
@@ -104,7 +104,7 @@ def plot_spectral_density(model, x_min=0.0, x_max=None):
 
 
 def plot_spectral_rad_pdf(model, x_min=0.0, x_max=None):
-    """plot radial spectral propability density function of a given CovModel"""
+    """Plot radial spectral pdf of a given CovModel."""
     if x_max is None:
         x_max = 3 / model.integral_scale
     x_s = np.linspace(x_min, x_max)
diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py
index 65c28b14e..c607a1d43 100644
--- a/gstools/covmodel/tools.py
+++ b/gstools/covmodel/tools.py
@@ -31,7 +31,7 @@
 else:
 
     class InitSubclassMeta(type):
-        """Metaclass that implements PEP 487 protocol
+        """Metaclass that implements PEP 487 protocol.
 
         Notes
         -----
@@ -43,6 +43,7 @@ class InitSubclassMeta(type):
         """
 
         def __new__(cls, name, bases, ns, **kwargs):
+            """Create a new subclass."""
             __init_subclass__ = ns.pop("__init_subclass__", None)
             if __init_subclass__:
                 __init_subclass__ = classmethod(__init_subclass__)
@@ -62,7 +63,7 @@ def __init__(cls, name, bases, ns, **kwargs):
 
 
 def rad_fac(dim, r):
-    """The volume element of the n-dimensional spherical coordinates.
+    """Volume element of the n-dimensional spherical coordinates.
 
     Given as a factor for integration of a radial-symmetric function.
 
@@ -90,7 +91,7 @@ def rad_fac(dim, r):
 
 
 def set_len_anis(dim, len_scale, anis):
-    """Setting the length scale and anisotropy factors for the given dimension.
+    """Set the length scale and anisotropy factors for the given dimension.
 
     Parameters
     ----------
@@ -139,7 +140,7 @@ def set_len_anis(dim, len_scale, anis):
                 constant_values=out_len_scale,
             )
         # if multiple length-scales are given, calculate the anisotropies
-        out_anis = np.zeros(dim - 1, dtype=float)
+        out_anis = np.zeros(dim - 1, dtype=np.double)
         for i in range(1, dim):
             out_anis[i - 1] = ls_tmp[i] / ls_tmp[0]
 
@@ -152,7 +153,7 @@ def set_len_anis(dim, len_scale, anis):
 
 
 def set_angles(dim, angles):
-    """Setting the angles for the given dimension.
+    """Set the angles for the given dimension.
 
     Parameters
     ----------
diff --git a/gstools/covmodel/tpl_models.py b/gstools/covmodel/tpl_models.py
index cc5a0a18d..72d6d6cce 100644
--- a/gstools/covmodel/tpl_models.py
+++ b/gstools/covmodel/tpl_models.py
@@ -26,7 +26,7 @@
 
 
 class TPLGaussian(CovModel):
-    r"""Truncated-Power-Law with Gaussian modes
+    r"""Truncated-Power-Law with Gaussian modes.
 
     Notes
     -----
@@ -108,14 +108,14 @@ class TPLGaussian(CovModel):
 
     @property
     def len_up(self):
-        """:class:`float`: The upper length scale truncation of the model.
+        """:class:`float`: Upper length scale truncation of the model.
 
             * ``len_up = len_low + len_scale``
         """
         return self.len_low + self.len_scale
 
     def var_factor(self):
-        r"""Factor for C (Power-Law factor) to result in variance
+        r"""Factor for C (Power-Law factor) to result in variance.
 
         This is used to result in the right variance, which is depending
         on the hurst coefficient and the length-scale extents
@@ -133,7 +133,7 @@ def var_factor(self):
         ) / (2 * self.hurst)
 
     def default_opt_arg(self):
-        """The defaults for the optional arguments:
+        """Defaults for the optional arguments.
 
             * ``{"hurst": 0.5, "len_low": 0.0}``
 
@@ -145,7 +145,7 @@ def default_opt_arg(self):
         return {"hurst": 0.5, "len_low": 0.0}
 
     def default_opt_arg_bounds(self):
-        """The defaults boundaries for the optional arguments:
+        """Defaults for boundaries of the optional arguments.
 
             * ``{"hurst": [0, 1, "oo"], "len_low": [0, 1000, "cc"]}``
 
@@ -157,7 +157,7 @@ def default_opt_arg_bounds(self):
         return {"hurst": [0, 1, "oo"], "len_low": [0, 1000, "cc"]}
 
     def correlation(self, r):
-        r"""Truncated-Power-Law with Gaussian modes - correlation function
+        r"""Truncated-Power-Law with Gaussian modes - correlation function.
 
         If ``len_low=0`` we have a simple representation:
 
@@ -196,7 +196,7 @@ def correlation(self, r):
 
 
 class TPLExponential(CovModel):
-    r"""Truncated-Power-Law with Exponential modes
+    r"""Truncated-Power-Law with Exponential modes.
 
     Notes
     -----
@@ -276,14 +276,14 @@ class TPLExponential(CovModel):
 
     @property
     def len_up(self):
-        """:class:`float`: The upper length scale truncation of the model.
+        """:class:`float`: Upper length scale truncation of the model.
 
             * ``len_up = len_low + len_scale``
         """
         return self.len_low + self.len_scale
 
     def var_factor(self):
-        r"""Factor for C (Power-Law factor) to result in variance
+        r"""Factor for C (Power-Law factor) to result in variance.
 
         This is used to result in the right variance, which is depending
         on the hurst coefficient and the length-scale extents
@@ -301,7 +301,7 @@ def var_factor(self):
         ) / (2 * self.hurst)
 
     def default_opt_arg(self):
-        """The defaults for the optional arguments:
+        """Defaults for the optional arguments.
 
             * ``{"hurst": 0.5, "len_low": 0.0}``
 
@@ -313,7 +313,7 @@ def default_opt_arg(self):
         return {"hurst": 0.5, "len_low": 0.0}
 
     def default_opt_arg_bounds(self):
-        """The defaults boundaries for the optional arguments:
+        """Defaults for boundaries of the optional arguments.
 
             * ``{"hurst": [0, 1, "oo"], "len_low": [0, 1000, "cc"]}``
 
@@ -325,7 +325,7 @@ def default_opt_arg_bounds(self):
         return {"hurst": [0, 1, "oo"], "len_low": [0, 1000, "cc"]}
 
     def correlation(self, r):
-        r"""Truncated-Power-Law with Exponential modes - correlation function
+        r"""Truncated-Power-Law with Exponential modes - correlation function.
 
         If ``len_low=0`` we have a simple representation:
 
@@ -362,7 +362,7 @@ def correlation(self, r):
 
 
 class TPLStable(CovModel):
-    r"""Truncated-Power-Law with Stable modes
+    r"""Truncated-Power-Law with Stable modes.
 
     Notes
     -----
@@ -453,14 +453,14 @@ class TPLStable(CovModel):
 
     @property
     def len_up(self):
-        """:class:`float`: The upper length scale truncation of the model.
+        """:class:`float`: Upper length scale truncation of the model.
 
             * ``len_up = len_low + len_scale``
         """
         return self.len_low + self.len_scale
 
     def var_factor(self):
-        r"""Factor for C (Power-Law factor) to result in variance
+        r"""Factor for C (Power-Law factor) to result in variance.
 
         This is used to result in the right variance, which is depending
         on the hurst coefficient and the length-scale extents
@@ -478,7 +478,7 @@ def var_factor(self):
         ) / (2 * self.hurst)
 
     def default_opt_arg(self):
-        """The defaults for the optional arguments:
+        """Defaults for the optional arguments.
 
             * ``{"hurst": 0.5, "alpha": 1.5, "len_low": 0.0}``
 
@@ -490,7 +490,7 @@ def default_opt_arg(self):
         return {"hurst": 0.5, "alpha": 1.5, "len_low": 0.0}
 
     def default_opt_arg_bounds(self):
-        """The defaults boundaries for the optional arguments:
+        """Defaults for boundaries of the optional arguments.
 
             * ``{"hurst": [0, 1, "oo"], "alpha": [0, 2, "oc"], "len_low": [0, 1000, "cc"]}``
 
@@ -506,7 +506,7 @@ def default_opt_arg_bounds(self):
         }
 
     def check_opt_arg(self):
-        """Checks for the optional arguments
+        """Check the optional arguments.
 
         Warns
         -----
@@ -521,7 +521,7 @@ def check_opt_arg(self):
             )
 
     def correlation(self, r):
-        r"""Truncated-Power-Law with Stable modes - correlation function
+        r"""Truncated-Power-Law with Stable modes - correlation function.
 
         If ``len_low=0`` we have a simple representation:
 
diff --git a/gstools/field/generator.py b/gstools/field/generator.py
index caa423e3d..4ec29e809 100644
--- a/gstools/field/generator.py
+++ b/gstools/field/generator.py
@@ -86,7 +86,7 @@ def __init__(
         self.update(model, seed)
 
     def __call__(self, x, y=None, z=None, mesh_type="unstructured"):
-        """Calculates the random modes for the randomization method.
+        """Calculate the random modes for the randomization method.
 
         This method  calls the `summate_*` Cython methods, which are the
         heart of the randomization method.
@@ -239,7 +239,7 @@ def reset_seed(self, seed=np.nan):
 
     @property
     def seed(self):
-        """:class:`int`: the seed of the master RNG
+        """:class:`int`: Seed of the master RNG.
 
         Notes
         -----
@@ -255,8 +255,7 @@ def seed(self, new_seed):
 
     @property
     def model(self):
-        """:any:`CovModel`: The covariance model of the spatial random field.
-        """
+        """:any:`CovModel`: Covariance model of the spatial random field."""
         return self._model
 
     @model.setter
@@ -265,7 +264,7 @@ def model(self, model):
 
     @property
     def mode_no(self):
-        """:class:`int`: The number of modes in the randomization method."""
+        """:class:`int`: Number of modes in the randomization method."""
         return self._mode_no
 
     @mode_no.setter
@@ -276,7 +275,7 @@ def mode_no(self, mode_no):
 
     @property
     def verbose(self):
-        """:class:`bool`: verbosity of the generator"""
+        """:class:`bool`: Verbosity of the generator."""
         return self._verbose
 
     @verbose.setter
@@ -285,20 +284,22 @@ def verbose(self, verbose):
 
     @property
     def name(self):
-        """:class:`str`: The name of the generator"""
+        """:class:`str`: Name of the generator."""
         return self.__class__.__name__
 
     def __str__(self):
+        """Return String representation."""
         return self.__repr__()
 
     def __repr__(self):
+        """Return String representation."""
         return "RandMeth(model={0}, mode_no={1}, seed={2})".format(
             repr(self.model), self._mode_no, self.seed
         )
 
 
 class IncomprRandMeth(RandMeth):
-    r"""Overrides RandMeth for calculating isotropic spatial incompressible random vector fields.
+    r"""RandMeth for incompressible random vector fields.
 
     Parameters
     ----------
@@ -338,8 +339,8 @@ class IncomprRandMeth(RandMeth):
         * :math:`Z_{k,j}` : random samples from a normal distribution
         * :math:`k_j` : samples from the spectral density distribution of
           the covariance model
-        * :math:`p_i(k_j) = e_1 - \frac{k_i k_1}{k^2}` : the projector ensuring the
-          incompressibility
+        * :math:`p_i(k_j) = e_1 - \frac{k_i k_1}{k^2}` : the projector
+          ensuring the incompressibility
     """
 
     def __init__(
@@ -353,7 +354,8 @@ def __init__(
     ):
         if model.dim < 2:
             raise ValueError(
-                "Only 2- and 3-dimensional incompressible fields can be generated."
+                "Only 2- and 3-dimensional incompressible fields "
+                + "can be generated."
             )
         super(IncomprRandMeth, self).__init__(
             model, mode_no, seed, verbose, **kwargs
@@ -362,10 +364,11 @@ def __init__(
         self.mean_u = mean_velocity
 
     def __call__(self, x, y=None, z=None, mesh_type="unstructured"):
-        """Overrides the Calculation of the random modes for the randomization method.
+        """Calculation of the random modes for the randomization method.
 
-        This method  calls the `summate_incompr_*` Cython methods, which are the
-        heart of the randomization method. In this class the method contains a projector to
+        This method  calls the `summate_incompr_*` Cython methods,
+        which are the heart of the randomization method.
+        In this class the method contains a projector to
         ensure the incompressibility of the vector field.
 
         Parameters
@@ -411,13 +414,15 @@ def __call__(self, x, y=None, z=None, mesh_type="unstructured"):
         )
 
     def _create_unit_vector(self, broadcast_shape, axis=0):
-        """Creates a unit vector which can be multiplied with a vector of shape broadcast_shape
+        """Create a unit vector.
+
+        Can be multiplied with a vector of shape broadcast_shape
 
         Parameters
         ----------
         broadcast_shape : :class:`tuple`
-            the shape of the array with which the unit vector is to be multiplied
-
+            the shape of the array with which
+            the unit vector is to be multiplied
         axis : :class:`int`, optional
             the direction of the unit vector. Default: ``0``
 
diff --git a/gstools/field/srf.py b/gstools/field/srf.py
index 123972786..f0e2014b7 100644
--- a/gstools/field/srf.py
+++ b/gstools/field/srf.py
@@ -235,28 +235,26 @@ def del_condition(self):
 
     @property
     def cond_pos(self):
-        """:class:`list`: The position tuple of the conditions.
-        """
+        """:class:`list`: The position tuple of the conditions."""
         return self._cond_pos
 
     @property
     def cond_val(self):
-        """:class:`list`: The values of the conditions.
-        """
+        """:class:`list`: The values of the conditions."""
         return self._cond_val
 
     @property
     def condition(self):
-        """:any:`bool`: State if conditions ar given"""
+        """:any:`bool`: State if conditions ar given."""
         return self._cond_pos is not None
 
     def cond_func(self, *args, **kwargs):
-        """The conditioning method applied to the field"""
+        """Conditioning method applied to the field."""
         if self.condition:
             return CONDITION[self._krige_type](*args, **kwargs)
 
     def structured(self, *args, **kwargs):
-        """Generate an SRF on a structured mesh
+        """Generate an SRF on a structured mesh.
 
         See :any:`SRF.__call__`
         """
@@ -264,7 +262,7 @@ def structured(self, *args, **kwargs):
         return call(*args, **kwargs)
 
     def unstructured(self, *args, **kwargs):
-        """Generate an SRF on an unstructured mesh
+        """Generate an SRF on an unstructured mesh.
 
         See :any:`SRF.__call__`
         """
@@ -272,11 +270,11 @@ def unstructured(self, *args, **kwargs):
         return call(*args, **kwargs)
 
     def upscaling_func(self, *args, **kwargs):
-        """The upscaling method applied to the field variance"""
+        """Upscaling method applied to the field variance."""
         return self._upscaling_func(*args, **kwargs)
 
     def set_generator(self, generator, **generator_kwargs):
-        """Set the generator for the field
+        """Set the generator for the field.
 
         Parameters
         ----------
@@ -302,8 +300,7 @@ def generator(self):
 
     @property
     def upscaling(self):  # pragma: no cover
-        """:class:`str`: Name of the upscaling method for the variance at each
-        point depending on the related element volume.
+        """:class:`str`: Name of the upscaling method.
 
         See the ``point_volumes`` keyword in the :any:`SRF.__call__` routine.
         Default: "no_scaling"
@@ -322,8 +319,7 @@ def upscaling(self, upscaling):
 
     @property
     def model(self):
-        """:any:`CovModel`: The covariance model of the spatial random field.
-        """
+        """:any:`CovModel`: The covariance model of the random field."""
         return self._model
 
     @model.setter
@@ -337,15 +333,15 @@ def model(self, model):
 
     @property
     def do_rotation(self):
-        """:any:`bool`: State if a rotation should be performed
-        depending on the model.
-        """
+        """:any:`bool`: State if a rotation is performed."""
         return not np.all(np.isclose(self.model.angles, 0.0))
 
     def __str__(self):
+        """Return String representation."""
         return self.__repr__()
 
     def __repr__(self):
+        """Return String representation."""
         return "SRF(model={0}, mean={1}, generator={2}".format(
             self.model, self.mean, self.generator
         )
diff --git a/gstools/field/upscaling.py b/gstools/field/upscaling.py
index 486d91cec..912f8b038 100644
--- a/gstools/field/upscaling.py
+++ b/gstools/field/upscaling.py
@@ -22,7 +22,7 @@
 
 
 def var_coarse_graining(model, point_volumes=0.0):
-    r"""Coarse Graning procedure to upscale the variance for uniform flow
+    r"""Coarse Graning procedure to upscale the variance for uniform flow.
 
     Parameters
     ----------
@@ -81,7 +81,7 @@ def var_coarse_graining(model, point_volumes=0.0):
 
 
 def var_no_scaling(model, *args, **kwargs):
-    r"""Dummy function to bypass scaling
+    r"""Dummy function to bypass scaling.
 
     Parameters
     ----------
diff --git a/gstools/random/rng.py b/gstools/random/rng.py
index db3c40b98..94dad1587 100644
--- a/gstools/random/rng.py
+++ b/gstools/random/rng.py
@@ -34,7 +34,6 @@ class RNG(object):
     """
 
     def __init__(self, seed=None):
-        """Initialize a random number generator"""
         # set seed
         self._master_rng = None
         self.seed = seed
@@ -48,7 +47,7 @@ def sample_ln_pdf(
         burn_in=20,
         oversampling_factor=10,
     ):
-        """Sample from a distribution given by ln(pdf)
+        """Sample from a distribution given by ln(pdf).
 
         This algorithm uses the :any:`emcee.EnsembleSampler`
 
@@ -127,7 +126,7 @@ def sample_ln_pdf(
         return self.random.choice(samples, size)
 
     def sample_dist(self, pdf=None, cdf=None, ppf=None, size=None, **kwargs):
-        """Sample from a distribution given by pdf, cdf and/or ppf
+        """Sample from a distribution given by pdf, cdf and/or ppf.
 
         Parameters
         ----------
@@ -163,7 +162,7 @@ def sample_dist(self, pdf=None, cdf=None, ppf=None, size=None, **kwargs):
         return dist.rvs(size=size)
 
     def sample_sphere(self, dim, size=None):
-        """Uniform sampling on a d-dimensional sphere
+        """Uniform sampling on a d-dimensional sphere.
 
         Parameters
         ----------
@@ -178,9 +177,9 @@ def sample_sphere(self, dim, size=None):
             x[, y[, z]] coordinates on the sphere with shape (dim, size)
         """
         if size is None:  # pragma: no cover
-            coord = np.empty(dim, dtype=float)
+            coord = np.empty(dim, dtype=np.double)
         else:
-            coord = np.empty((dim, size), dtype=float)
+            coord = np.empty((dim, size), dtype=np.double)
         if dim == 1:
             coord[0] = self.random.choice([-1, 1], size=size)
         elif dim == 2:
@@ -197,9 +196,9 @@ def sample_sphere(self, dim, size=None):
 
     @property
     def random(self):
-        """:any:`numpy.random.RandomState`:
-        Get a stream to the numpy Random number generator
+        """:any:`numpy.random.RandomState`: Randomstate.
 
+        Get a stream to the numpy Random number generator.
         You can use this, to call any provided distribution
         from :any:`numpy.random.RandomState`.
         """
@@ -207,7 +206,7 @@ def random(self):
 
     @property  # pragma: no cover
     def seed(self):
-        """:class:`int`: the seed of the master RNG
+        """:class:`int`: Seed of the master RNG.
 
         The setter property not only saves the new seed, but also creates
         a new master RNG function with the new seed.
@@ -219,9 +218,11 @@ def seed(self, new_seed=None):
         self._master_rng = MasterRNG(new_seed)
 
     def __str__(self):
+        """Return String representation."""
         return self.__repr__()
 
     def __repr__(self):
+        """Return String representation."""
         return "RNG(seed={})".format(self.seed)
 
 
diff --git a/gstools/random/tools.py b/gstools/random/tools.py
index 8e7f58e0b..361c1ab32 100644
--- a/gstools/random/tools.py
+++ b/gstools/random/tools.py
@@ -39,7 +39,7 @@ def __call__(self):
 
     @property  # pragma: no cover
     def seed(self):
-        """:class:`int`: the seed of the master RNG
+        """:class:`int`: Seed of the master RNG.
 
         The setter property not only saves the new seed, but also creates
         a new master RNG function with the new seed.
@@ -47,14 +47,16 @@ def seed(self):
         return self._seed
 
     def __str__(self):
+        """Return String representation."""
         return self.__repr__()
 
     def __repr__(self):
+        """Return String representation."""
         return "RNG(seed={})".format(self.seed)
 
 
 def dist_gen(pdf_in=None, cdf_in=None, ppf_in=None, **kwargs):
-    """Distribution Factory
+    """Distribution Factory.
 
     Parameters
     ----------
@@ -101,7 +103,7 @@ def dist_gen(pdf_in=None, cdf_in=None, ppf_in=None, **kwargs):
 
 
 class DistPdf(rv_continuous):
-    "Generate distribution from pdf"
+    """Generate distribution from pdf."""
 
     def __init__(self, pdf_in, **kwargs):
         self.pdf_in = pdf_in
@@ -112,7 +114,7 @@ def _pdf(self, x, *args):
 
 
 class DistCdf(rv_continuous):
-    "Generate distribution from cdf"
+    """Generate distribution from cdf."""
 
     def __init__(self, cdf_in, **kwargs):
         self.cdf_in = cdf_in
@@ -123,7 +125,7 @@ def _cdf(self, x, *args):
 
 
 class DistPdfCdf(rv_continuous):
-    "Generate distribution from pdf and cdf"
+    """Generate distribution from pdf and cdf."""
 
     def __init__(self, pdf_in, cdf_in, **kwargs):
         self.pdf_in = pdf_in
@@ -138,7 +140,7 @@ def _cdf(self, x, *args):
 
 
 class DistPdfPpf(rv_continuous):
-    "Generate distribution from pdf and ppf"
+    """Generate distribution from pdf and ppf."""
 
     def __init__(self, pdf_in, ppf_in, **kwargs):
         self.pdf_in = pdf_in
@@ -153,7 +155,7 @@ def _ppf(self, q, *args):
 
 
 class DistCdfPpf(rv_continuous):
-    "Generate distribution from cdf and ppf"
+    """Generate distribution from cdf and ppf."""
 
     def __init__(self, cdf_in, ppf_in, **kwargs):
         self.cdf_in = cdf_in
@@ -168,7 +170,7 @@ def _ppf(self, q, *args):
 
 
 class DistPdfCdfPpf(rv_continuous):
-    "Generate distribution from pdf, cdf and ppf"
+    """Generate distribution from pdf, cdf and ppf."""
 
     def __init__(self, pdf_in, cdf_in, ppf_in, **kwargs):
         self.pdf_in = pdf_in
diff --git a/gstools/tools/export.py b/gstools/tools/export.py
index 6aecc0833..de251689b 100644
--- a/gstools/tools/export.py
+++ b/gstools/tools/export.py
@@ -25,7 +25,7 @@
 
 
 def vtk_export_structured(filename, pos, field, fieldname="field"):
-    """Export a field to vtk structured rectilinear grid file
+    """Export a field to vtk structured rectilinear grid file.
 
     Parameters
     ----------
@@ -56,7 +56,7 @@ def vtk_export_structured(filename, pos, field, fieldname="field"):
 
 
 def vtk_export_unstructured(filename, pos, field, fieldname="field"):
-    """Export a field to vtk structured rectilinear grid file
+    """Export a field to vtk structured rectilinear grid file.
 
     Parameters
     ----------
@@ -88,7 +88,7 @@ def vtk_export_unstructured(filename, pos, field, fieldname="field"):
 def vtk_export(
     filename, pos, field, fieldname="field", mesh_type="unstructured"
 ):
-    """Export a field to vtk
+    """Export a field to vtk.
 
     Parameters
     ----------
diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py
index a2241be0f..0ec2116b8 100644
--- a/gstools/tools/geometric.py
+++ b/gstools/tools/geometric.py
@@ -82,7 +82,7 @@ def r3d_z(theta):
 
 
 def pos2xyz(pos, dtype=None, calc_dim=False):
-    """Convert postional arguments to x, y, z
+    """Convert postional arguments to x, y, z.
 
     Parameters
     ----------
@@ -126,7 +126,7 @@ def pos2xyz(pos, dtype=None, calc_dim=False):
 
 
 def xyz2pos(x, y=None, z=None, dtype=None):
-    """Convert postional arguments to x, y, z
+    """Convert postional arguments to x, y, z.
 
     Parameters
     ----------
diff --git a/gstools/tools/special.py b/gstools/tools/special.py
index 874c7836d..0c31c79b3 100644
--- a/gstools/tools/special.py
+++ b/gstools/tools/special.py
@@ -25,7 +25,7 @@
 
 
 def inc_gamma(s, x):
-    r"""The (upper) incomplete gamma function
+    r"""The (upper) incomplete gamma function.
 
     Given by: :math:`\Gamma(s,x) = \int_x^{\infty} t^{s-1}\,e^{-t}\,{\rm d}t`
 
@@ -46,7 +46,7 @@ def inc_gamma(s, x):
 
 
 def exp_int(s, x):
-    r"""The exponential integral :math:`E_s(x)`
+    r"""The exponential integral :math:`E_s(x)`.
 
     Given by: :math:`E_s(x) = \int_1^\infty \frac{e^{-xt}}{t^s}\,\mathrm dt`
 
@@ -65,7 +65,7 @@ def exp_int(s, x):
 
 
 def stable_cov_norm(r, len_scale, hurst, alpha):
-    r"""The normalized covariance function of the stable model
+    r"""The normalized covariance function of the stable model.
 
     Given by
 
@@ -87,7 +87,7 @@ def stable_cov_norm(r, len_scale, hurst, alpha):
     alpha : :class:`float`, optional
         Shape parameter of the stable model.
     """
-    r = np.array(np.abs(r / len_scale), dtype=float)
+    r = np.array(np.abs(r / len_scale), dtype=np.double)
     r[r < 1e-8] = 0  # hack to prevent numerical errors
     res = np.ones_like(r)
     res[r > 0] = (2 * hurst / alpha) * exp_int(
@@ -97,7 +97,7 @@ def stable_cov_norm(r, len_scale, hurst, alpha):
 
 
 def inc_beta(a, b, x):
-    r"""The incomplete Beta function
+    r"""The incomplete Beta function.
 
     Given by: :math:`B(a,b;\,x) = \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt`
 
diff --git a/gstools/tools/transform.py b/gstools/tools/transform.py
index 756dda6e6..db1138e8f 100644
--- a/gstools/tools/transform.py
+++ b/gstools/tools/transform.py
@@ -33,7 +33,7 @@
 
 def zinnharvey(field, conn="high", mean=None, var=None):
     """
-    Zinn and Harvey transformation to connect low or high values
+    Zinn and Harvey transformation to connect low or high values.
 
     Parameters
     ----------
@@ -97,7 +97,7 @@ def normal_force_moments(field, mean=0, var=1):
 
 def normal_to_lognormal(field):
     """
-    Transform normal distribution to log-normal distribution
+    Transform normal distribution to log-normal distribution.
 
     Parameters
     ----------
@@ -115,7 +115,7 @@ def normal_to_lognormal(field):
 
 def normal_to_uniform(field, mean=None, var=None):
     """
-    Transform normal distribution to uniform distribution on [0, 1]
+    Transform normal distribution to uniform distribution on [0, 1].
 
     Parameters
     ----------
@@ -144,7 +144,7 @@ def normal_to_uniform(field, mean=None, var=None):
 
 def normal_to_arcsin(field, mean=None, var=None, a=0, b=1):
     """
-    Transform normal distribution to arcsin distribution
+    Transform normal distribution to arcsin distribution.
 
     See: https://en.wikipedia.org/wiki/Arcsine_distribution
 
@@ -175,7 +175,7 @@ def normal_to_arcsin(field, mean=None, var=None, a=0, b=1):
 
 def normal_to_uquad(field, mean=None, var=None, a=0, b=1):
     """
-    Transform normal distribution to U-quadratic distribution
+    Transform normal distribution to U-quadratic distribution.
 
     See: https://en.wikipedia.org/wiki/U-quadratic_distribution
 
@@ -206,7 +206,7 @@ def normal_to_uquad(field, mean=None, var=None, a=0, b=1):
 
 def _uniform_to_arcsin(field, a=0, b=1):
     """
-    PPF of your desired distribution
+    PPF of your desired distribution.
 
     The PPF is the inverse of the CDF and is used to sample a distribution
     from uniform distributed values on [0, 1]
@@ -219,7 +219,7 @@ def _uniform_to_arcsin(field, a=0, b=1):
 
 def _uniform_to_uquad(field, a=0, b=1):
     """
-    PPF of your desired distribution
+    PPF of your desired distribution.
 
     The PPF is the inverse of the CDF and is used to sample a distribution
     from uniform distributed values on [0, 1]
diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py
index 92e60f16a..250414cce 100644
--- a/gstools/variogram/variogram.py
+++ b/gstools/variogram/variogram.py
@@ -122,7 +122,6 @@ def vario_estimate_structured(field, direction="x"):
     :class:`numpy.ndarray`
         the estimated variogram along the given direction.
     """
-
     try:
         mask = np.array(field.mask, dtype=np.int32)
         field = np.ma.array(field, ndmin=1, dtype=np.double)