code: def get_geostationary_bounding_box(geos_area, nb_points=50):
      
      xmax, ymax = get_geostationary_angle_extent(geos_area)
  
      # generate points around the north hemisphere in satellite projection
      # make it a bit smaller so that we stay inside the valid area
      x = np.cos(np.linspace(-np.pi, 0, nb_points / 2)) * (xmax - 0.001)
      y = -np.sin(np.linspace(-np.pi, 0, nb_points / 2)) * (ymax - 0.001)
  
      # clip the projection coordinates to fit the area extent of geos_area
      ll_x, ll_y, ur_x, ur_y = (np.array(geos_area.area_extent) /
                                geos_area.proj_dict['h'])
  
      x = np.clip(np.concatenate([x, x[::-1]]), min(ll_x, ur_x), max(ll_x, ur_x))
      y = np.clip(np.concatenate([y, -y]), min(ll_y, ur_y), max(ll_y, ur_y))
  
      return _lonlat_from_geos_angle(x, y, geos_area)
docs: Get the bbox in lon/lats of the valid pixels inside *geos_area*.
  
      Args:
        nb_points: Number of points on the polygon