Transforming an (lat,lon) extent in GOES-16 projection with Python/cartopy/GDAL?

Question

I have a WGS84 shapefile, which I'm trying to get its extent, firstly with its projection:

gdf = gpd.read_file("/content/buffer.shp")

bounds = row.geometry.bounds 

extent = [bounds[0], bounds[1], bounds[2], bounds[3]]
extent

Now I would like to search this extent in an GOES-16 .nc file. However, my shape is in WGS84 and GOES-16 is in a Geostationary projection. As I will search this extent in many .nc files, it would be easier to transform this shapefile extent rather than all .nc files to WGS84.

Can someone help me to do that?

I tried:

import cartopy.crs as ccrs
import xarray

def transform_extent(extent, source_proj, target_proj):
    transformer = pyproj.Transformer.from_crs(source_proj, target_proj, always_xy=True)
    min_lon, min_lat = transformer.transform(extent[0], extent[1])
    max_lon, max_lat = transformer.transform(extent[2], extent[3])
    return [min_lon, min_lat, max_lon, max_lat]

geo_proj = ccrs.Geostationary(central_longitude=-75.0, satellite_height=35786023.0).proj4_init
target_prj = osr.SpatialReference()
target_prj.ImportFromProj4(geo_proj)
transformed_extent = transform_extent(extent, 'EPSG:4326', geo_proj)

print(transformed_extent)


root='/content/GOES-16/ABI-L2-CMIPF/2019/001/07/OR_ABI-L2-CMIPF-M3C07_G16_s20190010700366_e20190010711145_c20190010711203.nc'


ds = xr.open_dataset(root)

perspective_point_height = file.variables['goes_imager_projection'].perspective_point_height
new_y = ds['y'][:] * perspective_point_height
new_x = ds['x'][:] * perspective_point_height

ds = ds.assign_coords(y=new_y, x=new_x)
min_lon, min_lat, max_lon, max_lat = x_min, y_min, x_max, y_max

min_lon=transformed_extent[0]
min_lat=transformed_extent[1]
max_lon=transformed_extent[2]
max_lat=transformed_extent[3]
cropped_ds = ds.sel(y=slice(max_lat, min_lat), x=slice(min_lon, max_lon))

That is, I tried to get y and x of the nc file (but I don't know if it is entirely correct by only multiplicating by perspective_point_height), tried to transform WGS84 -> Geoestationary and then do a slice. Is this workflow correct? Did I miss some projection step? Any easier way?

Any help is appreciated.