2

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.

2
0

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()

4
  • but only the cutline part is taking like 20 minutes, I am looking for a more effective way – Flash Thunder May 12 at 6:16
  • I tested on pretty large rasters with a small overlap (~300 x 10000 pixels overlap) and it takes seconds even accessing remote files via HTTP. Maybe the issue is that your TIFFs are not TILED? Even that is unlikely to be the case. Note that you mention you do not need to reproject here, so this is just data extraction. – Jose May 12 at 13:36
  • 48000 x 48000 px per stripe / about 450 000 x 450 000 px per Poland – Flash Thunder May 12 at 14:17
  • IO will take a while, but make sure your input data (img1 and img2) are tiled, and with a simple compression. At any case, my code should be reasonably fast as it is. – Jose May 12 at 15:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.