dlr-eoc / ukis-pysat

The ukis-pysat package provides generic classes and functions to query, access and process multi-spectral and SAR satellite images
Apache License 2.0
27 stars 8 forks source link

Bounds option for Image.warp #63

Closed MWieland closed 4 years ago

MWieland commented 4 years ago

It would be great to have an option to limit warping to a subset of the image in order to speed up computation. Currently, warping is always performed on the whole Image. An simple solution could be to add an optional "bounds" parameter that passes the bounding box coordinates in source crs coordinates to the method. A simple implementation could look as follows (NOTE: this example already considers #60 and #62 ):

    def warp(
        self, dst_crs, resampling_method=0, num_threads=4, resolution=None, bounds=None, nodata=None, target_align=None
    ):
        """Reproject a source raster to a destination raster.
        :param dst_crs: CRS or dict, Target coordinate reference system.
        :param resampling_method: Resampling algorithm, int, defaults to 0 (Nearest)
            numbers: https://github.com/mapbox/rasterio/blob/master/rasterio/enums.py#L28
        :param num_threads: int, number of workers, optional (default: 4)
        :param resolution: tuple (x resolution, y resolution) or float, optional.
            Target resolution, in units of target coordinate reference system.
        :param bounds: tuple bounding box coordinates in source coordinate reference system, optional.
        :param target_align: raster to which to align resolution, extent and gridspacing, optional (Image).
        :param nodata: nodata value of source, int or float, optional.
        """
        if target_align:
            transform = target_align.dataset.transform
            width = target_align.dataset.width
            height = target_align.dataset.height

        else:
            if not bounds:
                bounds = self.dataset.bounds

            if resolution:
                transform, width, height = calculate_default_transform(
                    self.dataset.crs, dst_crs, self.dataset.width, self.dataset.height, *bounds, resolution=resolution,
                )

            else:
                transform, width, height = calculate_default_transform(
                    self.dataset.crs, dst_crs, self.dataset.width, self.dataset.height, *bounds,
                )

        destination = np.zeros((self.dataset.count, height, width), self.__arr.dtype)

        self.__arr, transform = reproject(
            source=self.__arr,
            destination=destination,
            src_transform=self.dataset.transform,
            src_crs=self.dataset.crs,
            src_nodata=nodata,
            dst_transform=transform,
            dst_crs=dst_crs,
            dst_nodata=nodata,
            resampling=resampling_method,
            num_threads=num_threads,
        )

        self.__update_dataset(dst_crs, transform, nodata=nodata)
fwfichtner commented 4 years ago

Please be careful with the order of parameters to keep compatibility. I am also not sure how useful this is, why would you read in the whole image and then only use a subset of it? Would it not be better to read in a subset from the start?

MWieland commented 4 years ago

To work on a subset of an image can be very useful and efficient. I agree with you that reading in a subset from the start would be the better option. @fwfichtner Can we add a subset option (like windowed reading) to Image initialization?

fwfichtner commented 4 years ago

It is quite simple and already there if you want. You can read in the subset like in this example and the initialize Image with the windowed dataset.

MWieland commented 4 years ago

Thanks for the example. It seems to work with image coordinates though. Wouldnt it be desirable to have this work with real world coordinates? We could provide this as a convenience call option in ukis-pysat, which could proof an added value compared to plain rasterio windowed reading. Curious about your opinion.

fwfichtner commented 4 years ago

This you can solve using rasterio.windows.from_bounds. I guess what I'm saying is that we would just put something inside which already exists outside and is more accessible there.

MWieland commented 4 years ago

Alright in this case I agree. Did not know about this rasterio method. Thanks for clarifying. Closing this issue.