Python fastest way to load overlapping part of 2 rasters with different origin to np.array

Question

I got a problem with gdal.wrap - it works REALLY slow. On 10m raster of Poland, matching the origin of 2 rasters takes like 30 minutes.

Is there an easy and fast way to take 2 rasters (of different origins) and to load their overlapping part into the numpy arrays?

Rasters are having the same resolution and crs.

    def run_mosaic_script(self, resample_alg: any = None, x_res=13.9, y_res=13.9, compress_algorithm="LZW",
                          dstSRS="EPSG:32634") -> bool:
        """
        Method for run mosaic script
​
        :param compress_algorithm: compression algorithm example LZW
        :type compress_algorithm: str or None
        :param resample_alg: resample algorithm
        :type resample_alg: any
        :param x_res: x resolution of pixel
        :type x_res: numeric (int or float) or None
        :param y_res: resolution of pixel
        :type y_res:numeric (int or float) or None
        :return: True if correct, False if incorrect
        :rtype: bool
        """
        creation_option: List[str] = ["BIGTIFF=YES"]
        if compress_algorithm is not None:
            creation_option.append("COMPRESS=" + compress_algorithm)
        gdal_warp_options = gdal.WarpOptions(srcNodata=0, dstNodata=0, multithread=True,
                                             warpOptions="NUM_THREADS=ALL_CPUS", xRes=x_res, yRes=y_res,
                                             dstSRS=dstSRS,
                                             resampleAlg=resample_alg, creationOptions=creation_option,
                                             cutlineDSName=self.shp_file.__str__(), cropToCutline=True,
                                             overviewLevel='AUTO')
        image_to_mosaic: list = list()
        for file_name in self.input_image_list:
            image_to_mosaic.append(file_name.__str__())
        gdal.Warp(destNameOrDestDS=self.output_image.__str__(),
                  srcDSOrSrcDSTab=image_to_mosaic,
                  options=gdal_warp_options)
        return True

I am just cropping it to some simple extent. Both of the rasters.

Answer 1

Basically, you want to:

1. query the rasters & get the geometries for them,
2. find the intersection,
3. use gdal.Warp with a cutline to crop the intersected bit.

The first step can be done easily with gdal.Info, as it basically returns the extent which can then be used to calculate the intersection. Once you have this, you have to export it to a GeoJSON file, and do the warping.

I was hoping to use in-memory vector datasets for point 2 above, but I can't figure out how to use them, which would save you using a temporary GeoJSON file.

Note that even if the two datasets are in different projections, you can still do this and use the last gdal.Warp calls to reproject to a common projection/pixel size.

import json
from osgeo import gdal, ogr, osr

def get_overlap(img1, img2):
    extent = gdal.Info(img1, format="json")['wgs84Extent']
    poly1 = ogr.CreateGeometryFromJson(json.dumps(extent))
    extent = gdal.Info(img2, format="json")['wgs84Extent']
    poly2 = ogr.CreateGeometryFromJson(json.dumps(extent))
    intersection = poly1.Intersection(poly2)
    gg = gdal.OpenEx(intersection.ExportToJson())
    ds = gdal.VectorTranslate('/tmp/output.json', srcDS=gg, 
                              format = 'GeoJSON', 
                              layerCreationOptions = ['RFC7946=YES',
                                                      'WRITE_BBOX=YES'
                                                     ]
                             )
    ds = None # Flush to disk
    # Returns a GeoJSOn with the overlap area as polygon.
    # Not tested what happens if files don't overlap!
    return "/tmp/output.json"
    
#img1 = blah blah
#img2 = blah blah

overlap = get_overlap(img1, img2)

img1_crop = gdal.Warp("", img1, format="MEM", cropToCutline=True,
                    cutlineDSName=overlap).ReadAsArray()
img2_crop = gdal.Warp("", img2, format="MEM", cropToCutline=True,
                    cutlineDSName=overlap).ReadAsArray()