docs: def lonlat_to_laea(lon, lat, lon0, lat0, f_e=0.0, f_n=0.0):
      """
      Converts vectors of longitude and latitude into Lambert Azimuthal
      Equal Area projection (km), with respect to an origin point
  
      :param numpy.ndarray lon:
          Longitudes
      :param numpy.ndarray lat:
          Latitude
      :param float lon0:
          Central longitude
      :param float lat0:
          Central latitude
      :param float f_e:
          False easting (km)
      :param float f_e:
          False northing (km)
      :returns:
          * easting (km)
          * northing (km)
      """
      lon = np.radians(lon)
      lat = np.radians(lat)
      lon0 = np.radians(lon0)
      lat0 = np.radians(lat0)
      q_0 = TO_Q(lat0)
      q_p = TO_Q(np.pi / 2.)
      q_val = TO_Q(lat)
      beta = np.arcsin(q_val / q_p)
      beta0 = np.arcsin(q_0 / q_p)
      r_q = WGS84["a"] * np.sqrt(q_p / 2.)
      dval = WGS84["a"] * (
          np.cos(lat0) / np.sqrt(1.0 - (WGS84["e2"] * (np.sin(lat0) ** 2.))) /
          (r_q * np.cos(beta0)))
      bval = r_q * np.sqrt(
          2. / (1.0 + (np.sin(beta0) * np.sin(beta)) + (np.cos(beta) *
                                                        np.cos(beta0) * np.cos(lon - lon0))))
      easting = f_e + ((bval * dval) * (np.cos(beta) * np.sin(lon - lon0)))
      northing = f_n + (bval / dval) * ((np.cos(beta0) * np.sin(beta)) -
                                        (np.sin(beta0) * np.cos(beta) * np.cos(lon - lon0)))
      return easting, northing