geotools.py

from osgeo import gdal
from osgeo import ogr
from pyproj import Transformer as T

def get_query_bbox(shp, query_srs=4326, return_all=False):
    ds = ogr.Open(shp)
    ly = ds.GetLayer(0)
    xt = ly.GetExtent()
    shp_srs = int(ly.GetSpatialRef().GetAuthorityCode(None))

    i_bbox = [xt[0], xt[2], xt[1], xt[3]]
    if shp_srs == query_srs:
        bbox = i_bbox
    else:
        tr = T.from_crs(shp_srs, query_srs, always_xy=True)
        lon_min, lat_min = tr.transform(xt[0], xt[2])
        lon_max, lat_max = tr.transform(xt[1], xt[3])
        bbox = [lon_min, lat_min, lon_max, lat_max]

    if return_all:
        return bbox, i_bbox, shp_srs
    else:
        return bbox