Calculating Intersection Extent of Two GeoTIFF Files in Python

Python and gdal, get intersect extend of two geotiff file

To get the intersecting extent of two GeoTIFF files using Python and the GDAL library, you can follow these steps:

1. Open the GeoTIFF files.
2. Retrieve the bounding boxes (extents) of both files.
3. Calculate the intersection of these bounding boxes.

Here's a step-by-step guide with code:

1. Install GDAL: If you haven't already, you need to install the GDAL library. You can do this using pip:

   pip install gdal
   

2. Python Code:

   from osgeo import gdal
   
   def get_extent(ds):
       """Get the extent of a GDAL dataset."""
       gt = ds.GetGeoTransform()
       minx = gt[0]
       maxx = gt[0] + (ds.RasterXSize * gt[1])
       miny = gt[3] + (ds.RasterYSize * gt[5])
       maxy = gt[3]
       return (minx, miny, maxx, maxy)
   
   def intersect_extent(extent1, extent2):
       """Calculate the intersection of two extents."""
       minx = max(extent1[0], extent2[0])
       miny = max(extent1[1], extent2[1])
       maxx = min(extent1[2], extent2[2])
       maxy = min(extent1[3], extent2[3])
       
       if minx < maxx and miny < maxy:
           return (minx, miny, maxx, maxy)
       else:
           return None  # No intersection
   
   # Open the GeoTIFF files
   ds1 = gdal.Open('path_to_first_geotiff.tif')
   ds2 = gdal.Open('path_to_second_geotiff.tif')
   
   # Get extents
   extent1 = get_extent(ds1)
   extent2 = get_extent(ds2)
   
   # Calculate intersection
   intersection = intersect_extent(extent1, extent2)
   
   if intersection:
       print(f"Intersection extent: {intersection}")
   else:
       print("No intersection found.")
   

Explanation:

1. get_extent(ds): This function takes a GDAL dataset and returns its extent as a tuple (minx, miny, maxx, maxy).

2. intersect_extent(extent1, extent2): This function takes two extents and calculates their intersection. If the extents do not intersect, it returns None.

3. Opening GeoTIFF files: The gdal.Open function is used to open the GeoTIFF files.

4. Getting extents: The get_extent function is called for both datasets to get their extents.

5. Calculating intersection: The intersect_extent function is used to find the intersection of the two extents.

6. Output: The intersection extent is printed if it exists; otherwise, a message indicating no intersection is printed.

Make sure to replace 'path_to_first_geotiff.tif' and 'path_to_second_geotiff.tif' with the actual paths to your GeoTIFF files.