geotools.py

from osgeo import gdal
from osgeo import ogr
from pyproj import Transformer as T
from shapely.geometry import box
from shapely.ops import transform
import rasterio
import geopandas

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
    
def check_extent_overlap(im,vect):
    im_ds = rasterio.open(im)
    im_srs = im_ds.crs
    im_extent = box(*im_ds.bounds)

    vect_ds = geopandas.read_file(vect)
    vect_srs = vect_ds.crs
    vect_extent = box(*vect_ds.geometry.total_bounds)

    if (im_srs != vect_srs):
        tr = T.from_crs(im_srs, vect_srs, always_xy=True)
        im_extent = transform(tr.transform, im_extent)
        
    return im_extent.intersects(vect_extent)