docs: def get_distance_matrix(self):
          """
          Compute and return distances between each pairs of points in the mesh.
  
          This method requires that the coordinate arrays are one-dimensional.
          NB: the depth of the points is ignored
  
          .. warning::
              Because of its quadratic space and time complexity this method
              is safe to use for meshes of up to several thousand points. For
              mesh of 10k points it needs ~800 Mb for just the resulting matrix
              and four times that much for intermediate storage.
  
          :returns:
              Two-dimensional numpy array, square matrix of distances. The matrix
              has zeros on main diagonal and positive distances in kilometers
              on all other cells. That is, value in cell (3, 5) is the distance
              between mesh's points 3 and 5 in km, and it is equal to value
              in cell (5, 3).
  
          Uses :func:`openquake.hazardlib.geo.geodetic.geodetic_distance`.
          """
          assert self.lons.ndim == 1
          distances = geodetic.geodetic_distance(
              self.lons.reshape(self.lons.shape + (1, )),
              self.lats.reshape(self.lats.shape + (1, )),
              self.lons,
              self.lats)
          return numpy.matrix(distances, copy=False)