code: def _grid_sfc_area(lon, lat, lon_bounds=None, lat_bounds=None):
      
      # Compute the bounds if not given.
      if lon_bounds is None:
          lon_bounds = _bounds_from_array(
              lon, internal_names.LON_STR, internal_names.LON_BOUNDS_STR)
      if lat_bounds is None:
          lat_bounds = _bounds_from_array(
              lat, internal_names.LAT_STR, internal_names.LAT_BOUNDS_STR)
      # Compute the surface area.
      dlon = _diff_bounds(utils.vertcoord.to_radians(lon_bounds, is_delta=True),
                          lon)
      sinlat_bounds = np.sin(utils.vertcoord.to_radians(lat_bounds,
                                                        is_delta=True))
      dsinlat = np.abs(_diff_bounds(sinlat_bounds, lat))
      sfc_area = dlon*dsinlat*(RADIUS_EARTH**2)
      # Rename the coordinates such that they match the actual lat / lon.
      try:
          sfc_area = sfc_area.rename(
              {internal_names.LAT_BOUNDS_STR: internal_names.LAT_STR,
               internal_names.LON_BOUNDS_STR: internal_names.LON_STR})
      except ValueError:
          pass
      # Clean up: correct names and dimension order.
      sfc_area = sfc_area.rename(internal_names.SFC_AREA_STR)
      sfc_area[internal_names.LAT_STR] = lat
      sfc_area[internal_names.LON_STR] = lon
      return sfc_area.transpose()
docs: Calculate surface area of each grid cell in a lon-lat grid.